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DOCUMENTATION OF THE BENSON DIESEL ENGINE SIMULATION PROGRAM 


Jon Van Gerpen* 

Propulsion Directorate 

U.S. Army Aviation Research and Technology Activity - AVSCOM 

Lewis Research Center 
Cleveland, Ohio 44135 


1.0 INTRODUCTION 


This report documents the Benson Diesel Engine Simulation Program and 
explains how it can be used to predict the performance of diesel engines. 

This program was obtained from the Garrett Turbine Engine Company but has been 
extensively modified since that time. The program is a thermodynamic 
simulation of the diesel engine cycle which uses a single zone combustion 
model. It can be used to predict the effect of changes in engine design and 
operating parameters such as valve timing, speed and boost pressure. The most 
significant change made to this program is the addition of a more detailed 
heat transfer model to predict metal part temperatures. 

This report contains a description of the sub-models used in the Benson 
program, a description of the input parameters and sample program runs. 


*Iowa State University, Mechanical Engineering Dept., Ames, Iowa 50010 
and Summer Faculty Fellow at Lewis Research Center. 
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2.0 OVERVIEW OF PROGRAM 

As shown in Figure 2.1, the Benson program uses a large number of 
subroutines to accomplish specialized tasks during the calculation procedure. 
However, most of the cycle calculations are performed in two subroutines: 
CYLNDR, which calculates the gas exchange processes and POWER which takes care 
of compression, combustion and expansion. The development of the energy 
balance and a discussion of how it is evaluated is given below. 


2. 1 Energy Balance - Benson Model 

The chemical equation for the combustion process in the diesel engine can 

be simplified to the following expression. 

C H + L n +m/2 J [o 2 + 3.77 N 2 ] 
n m <p 

. m TT „ , fTn + m/2] m]^ . 0 ^ fn + m/2] RI 

-> n 00 2 + £ H 2 0 + l L 1 - n - ^J0 2 + 3.77 1 L N 2 

where n and m are stoichiometric coefficients and <f> is the equivalence ratio. 

The first law of thermodynamics must be satisfied at all times as this 

reaction occurs in the engine. For the period when the intake and exhaust 

ports are closed, the first law can be written for an interval of time from t 

to t+dt, 

E 2 - E t = 6Q - 6W 

where E 2 and E t are given by 

E 2 = aj^ e (T 2 ) + a Q e Q (T 2 ) + a oo e oo ( T 2 ) + ^ 2 O e H 2 0^ T2 ^ + a C H e C H 2 ^ 

2 2 22 2 z zz n m n m 

El = + b 0 2 e 0 2 f T ‘) + b C» 2 e C0 2 ^ Tl ^ + + b C H ®C H 


n m n m 


The coefficients, and are the number of moles of species i present in 


the cylinder and 


e.(T) = e.(0) + [e.(T) - e^O)] 
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Figure 2.1 Flow Chart for Benson Program 
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where e^fO) is the partial molar internal energy of species i at absolute zero 
and [e^(T) - e^(0)] is the sensible internal energy difference between 
absolute zero and T. 


For convenience, the sum of the sensible terms can be defined as E 

E*(T) = ^ 2 [' N2 (T) - <^(0)] + a o 2 l>o 2 ( T > - 

+ WW 0 " e C0 2 (0)] + a H 2 o[ e H 2 0< T > - e H 2 0 (0 » 

Thus E 2 - Ej can be written as 

E 2 - E, = (E* - E*) * (a^-b^Je^fO) * (a^-b^Je^fO) 

+ (a C0 2 “ b C0 2 )e C0, (0) + (a H 2 0' b H,0> e H 2 0'°> 


+ (a C H " b C H > e C H <°> 
n m n m n m 

If dMj. moles of fuel are burned as the reaction progresses from t to t+dt, 

atom balances can be used to fix relations between the a. and b. coefficients. 

1 1 

a W 2 = b N 2 

a o 2 = V [n + 5 

a co 2 = b C0 2 + n 
a H 2 0 = b H 2 0 + 2 


a C H = b c H ” 
n m n m 


The constant volume heat of combustion with the water in the products 


present as vapor can be defined as 


Q = E ,(T ) - E (T ) 
v prod v s' reac v s' 


where E prod - n e^O^) + 2 e H 2 0^ T s^ 

E reac = [ n + m/2 ] e 0 2 ^ T s^ + ®C H 

^ n m 

T = reference temperature where Q is defined. 
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By combining the definition of the constant volume heat of combustion and the 
definition of the constant volume heat of combustion at absolute zero 


we can obtain 


AE = E JO) - E (0) 
o prod v 1 reac v ' 


% - 1E o * "l>C0 2 < T s> - e C0 2 (°>] * fC%0< T s> - 
- (” * |)[ e 0 2 < T S ) - 'o 2 ‘° )] - [®C H <V - e C H <°)3 

n m n m 

Substituting this into the equation for (Es-Ej) gives 

E 2 - E t = E* - eT - Q v dM f - (a O02 -b CO2 )[e CO2 (T s ) - *^(0)] 


^ a H 2 0~ b H 2 0^ e H 2 0^s^ “ e H 2 0^°^ 
+ ^ a 0 2 “ b 0 2 ^ e H 2 0^ T s^ " e H 2 0^°^ 

+ a CH ^ e C H {T2)_e C H (VI 


n m n m 


n m 


■ b C H [e C H (T ‘ )_e C H < T S )] 
n m n m n m 

Since the fuel is such a small fraction of the total mass in the 

cylinder, the fuel internal energy terms can be neglected and the terms 

evaluated at T can be written as 
s 

E 2 - E t = [E* - E 2 (T s )] - [E? - E?(T s )] - Q v dM f 
where E 1 (T s ) designates the sum of the sensible internal energy terms using 
the composition at state 1 and evaluating the partial molar internal energies 
at T . 


Finally, the first law of thermodynamics can be written as 

[E* - E 2 (T s )] - [E? - E?(T s )] - Q y dM f = 5Q - 5W 

for the period where the intake and exhaust ports are closed and when the 

ports are open the first law becomes 

[E* - E*(T )] - [E? - E?(T )] - Q dM. = 6Q - 6W + m h - m h „ 

L 2 zv s' J L 1 1V s' J v f in in out out 
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Each of the terms in this equation will be discussed in turn. 

2.2 Change in Internal Energy 

The internal energy is calculated using fourth order polynomial curve 
fits of the JANAF tables [1]. These curve fits provide the sensible enthalpy 
change between absolute zero and some temperature T. The values of internal 
energy are obtained from e = h - RT . 

2.3 Combustion 

The rate at which fuel is burned, and its energy released, is determined 
by the choice of combustion model. Currently, four combustion models may be 
selected for use in the Benson Program along with one user supplied model. 

The choice of model is determined by the values of the input parameters, NPCWF 
and DEL3. The values of these parameters corresponding to each model are 
given below. 


NPCWF = 0 and DEL3 > 0. 

NPCWF = 0 and DEL3 = 0. 

NPCWF = -1 
NPCWF = 1 
NPCWF = 2 

Each of the combustion models 


Wiebe combustion model 
AVL combustion model 
Whi tehouse-Way combustion model 
Watson combustion model 
User supplied model 
11 be discussed in turn. 


2.3.1 Wiebe Model 

The Wiebe model is a simple functional model consisting of a two 
parameter exponential function which has the flexibility to assume the shape 
of known burning rate curves in supercharged diesel engines. The following 
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expression is used to calculate the fraction of fuel burned. 


% Fuel Burned = 100. x [1 - exp(-6j t 2 )] 

where 6 t and 6 2 are parameters which determine the shape of the function 

CA - ACB 


r is a dimensionless time parameter = 
CA is the current crank angle 


HRD 


ACB is the crank angle at which combustion begins 
HRD is the combustion duration, equal to (ACF - ACB) 


2.3.2 The AVL Model 

The AVL model is the same as the Wiebe model except that AVL has 
recommended a way to specify the combustion duration. AVL suggests that the 
combustion duration parameter be chosen so that it is a linear function of the 
equivalence ratio. 

HRD = 15 + 115 (<#.) 

where if> is the fuel-air equivalence ratio. 


2.3.3 The Whitehouse-Way Model 

Of the models supplied with the Benson Program, the Whitehouse-Way [2,3] 
model is the only one that attempts to model the physics of the combustion 
event. The fuel injection rate is supplied as am input and the model relates 
the burning rate to this injection rate. The injection rate is specified by 
two parameters: the initial injection rate and the timing at which the maximum 
injection rate occurs. The injection rate is defined by the following 
function. 


7 



. TaA 2 - AAN(l) 
Injection Rate = AAN(l) + | — AANm X - L 

AA2 


AAN(2) 
(r - 1) 


" AAN(2) - 1 

where AAN(l) and AAN(l) are user specified injection rate parameters 
AA2 = 2 - IaAN(2) *AAN( 1 ) 


0 < t < AAN(2) 
AAN(2) < t < 1 


t = defined same as above under Wiebe 

Fuel is assumed to be prepared for combustion according to a relation of the 
form 


P = K M, 1 X M X P n m 
i u 0 2 

where x and m are empirical parameters which should be constant for a given 
problem. 

K is a parameter which needs to be determined for each individual 
problem. 

P is the fuel preparation rate. 

is the mass of fuel injected. 

M is the mass of fuel unburnt, 
u 

The burning rate is calculated from an exponential expression of the form: 

R = Sr — 2 [(P - R) dx • exp(-^|h 


where R is the burning rate 

"act" is an effective activation energy which must be specified 
P is the partial pressure of oxygen 

o 2 

N is the engine speed, Rev/min 

The integral is evaluated in the Benson program as a summation 


2(P-R)ACA 


and R is calculated using the value of the summation calculated up to the 
beginning of the time step. 
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Although the authors of this model have attempted to incorporate 
combustion theory into this model, it should be recognized that it is based on 
an understanding of diesel combustion that assumes droplet effects are 
important. Current understanding of diesel combustion assumes that the fuel 
vaporizes very soon after it is injected and that the combustion has more in 
common with a turbulent gaseous diffusion flame than droplet combustion. 


2.3.4 The Watson Model 

The Watson model [4] is a functional model, like the Wiebe model, which 
uses simple algebraic functions chosen for their ability to produce burning 
curves similar to those observed in diesel engines. The Watson model was 
developed to reproduce the ’’premixed spike” at the start of typical diesel 
heat release curve. Instead of a single exponential curve like the Wiebe 
model, the Watson model uses two curves which are superimposed. The first 
curve has a very high rate of combustion and simulates the premixed 
combustion. This curve dies out very quickly after the start of combustion 
and the second curve, which is an exponential exactly like the Wiebe model, 
becomes dominant. A weighting factor is used to determine how much emphasis 
is given to each of the curves. The Watson model is as follows. 


% Fuel Burned = 100* 


0 C u. 

(l - [1 - T Pl ] Pz ] + (1 - /3)[l - exp[-C di T da ] 


6 2 “ 5 3 

where j3 = weighting factor = 1. - 6^ <f> actr 

C = 2. + 1.25 X 10" 8 (actr • Nr • 60. ) 2 ' 4 
Pi 


C = 5000. 
P 2 

C, = 14.2 <t> 
di 


-0.644 
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0.25 


C, = 0.79 • C 
d 2 < 

actr = ignition delay period, which may be estimated from 

actr = 3.45 exp(2100./T ) P _102 

nr m 

where T and P are the mean temperature and pressure 
m m 

during the ignition delay period. 

$ = equivalence ratio 
Nr = engine speed, rev/sec 

6 , fig* ^3 = parameters chosen to vary shape of the model. 

(Watson suggests 5^ = 0.926, 6 ^ ~ 0.37 and 6 ^ = -0.26) 

The Watson model, because it can simulate premixed combustion, is 
appropriate for conditions where premixed combustion is important such as 
light loads and cold starting with naturally aspirated engines. This model is 
also appropriate when running on low cetane fuels. Under highly supercharged 
conditions the amount of premixed combustion will be small ( 0 = 0 ) so the 
Watson model reduces to the Wiebe model. 

2.3.5 User supplied model 

The user may provide a mass burning rate function of any arbitrary shape 
through the subroutine BURNRT. The subroutine should be written so that at a 
given value of crank position, designated by ANGLA, the amount of fuel burned 
in a one degree time step is returned as DWFUEL. 

2.4 Heat Transfer 

Two correlations are available for use in calculating the heat transfer 
in the Benson Program: The Annand [5] correlation and the Woschni [ 6 ] 
correlation. The Annand correlation takes the following form. 
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q/A = (Re) b (T - T w ) + c (T 4 - T^) 
where q = heat transfer 

A = Area of surface exposed to heat transfer 
k = Thermal conductivity 

p U D 

Re = Reynolds number = — — 

T = Gas Temperature 

T^ = Wall Temperature 

U = Mean Piston speed 

P 

p = gas density 

p = gas viscosity 

D = Bore of engine 

a,b,c = User supplied constants, called ANNA, ANNB and ANNC in the Benson 
Program. 

The constant, a, depends on the intensity of gas motion in the cylinder 

and Annand recommends a value between 0.35 and 0.8, with higher values 

corresponding to more intense motion. The constant, b, should be 0.7 and the 

constant c should be- 

during compression, c = 0. 

during combustion and expansion, 

-11 vw 

c = 3.266 X 10 — ~ j for diesel engine 

m - K 
-12 kW 

c = 4.286 X 10 j for spark ignition engine 

in - k 

Woschni’s correlation is of the following form. 

q/A = a | (Re) b (T - T w ) 
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where Re = Reynolds Number = ^ ^ — 
w = In-cylinder gas velocity 


= 6.18 U during scavenging period, m/s. 

P 

= 2.28 U during compression, m/s 
P 

V T 1 

= 2.28 U + 3.24 X 10 s y (p - p ) during combustion 
P Pi M P 


1 1 


and expansion, m/s 


V g = displacement volume. 

Pj. Vj, = pressure, volume and temperature at some point when the 
cylinder is closed for use in estimating the trapped mass. 

p = ins tantaneous cylinder pressure 

p = cylinder pressure for an unfired cycle at the same crank position 
and under the same conditions as p. 
a.b = User supplied constants, called ANNA and ANNB in the Benson 

program. Woschni recommends a value of 0.035 for constant a and 
0.8 for constant b. 

The principal difference between these two correlations is that Annand’s 
correlation separates the convective and radiative terms while Woschni’s 
correlation combines the convection and radiation in a single term. The user 
supplied parameters ANNA, ANNB and ANNC are used to select which correlation 
will be used in the program and to adjust the significance of heat transfer in 
the energy balance. When ANNC is equal to 0.0, the Woschni correlation will 
be used, and when ANNC is equal to a non-zero number, the Annand correlation 
will be used. ANNB is normally set to the respective author’s suggested 
value. ANNA is also normally set to the author's values unless data is 
available which indicate the total amount of heat loss for the cycle in which 
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case ANNA can be adjusted to obtain this amount. 


2.5 Work 

The work term is calculated from time step to time step by taking the 
average pressure for the time step multiplied by the volume change for the 
time step. 


6W = . (V 2 - V,) 


2.6 Mass Flow Terms 

The mass flow terms are calculated using subroutine MASSFL which 
calculates the mass entering or leaving the cylinder through the intake or 
exhaust ports. The techniques involved in this calculation are explained 
later in the discussion of that subroutine. After calculating the mass flow 
rates, the enthalpy is calculated using the polynomial curve-fits to the JANAF 
tables . 


13 



3.0 STEADY STATE HEAT TRANSFER ANALYSIS 


3. 1 Introduction 

The original version of the Benson Program received from Garrett contains 
a simple heat transfer analysis which calculates convective heat transfer 
coefficients using either the Annand or the Woschni correlations as described 
in section 2.4. The analysis described in this section has been added to the 
Benson program to allow prediction of the effect of varying engine design 
parameters and operating conditions on metal part temperatures. The analysis 
provides an estimate of the gas and coolant side wall temperatures for the 
piston crown, the cylinder wall, the head and the exhaust valve along with the 
top compression ring temperature. 

The model is one -dimensional and assumes the temperatures of the various 
components do not change during the cycle. The model has been developed for a 
single cylinder engine and can not distinguish between end cylinders or 
internal cylinders of a multi-cylinder engine. Submodels for the heat 
generated by piston/ring friction and for crevice volume heat transfer are 
included. 

The analysis requires that a number of critical engine dimensions be 
specified along with heat transfer path lengths. Although these heat transfer 
path lengths can be estimated, a high level of accuracy requires that 
experimental measurements be used to calculate these path lengths and thus 
"tune" the heat transfer analysis. Without experimental data, the absolute 
values of the temperatures provided by this analysis can be no more accurate 
than the path lengths supplied as input. However, even though the numerical 
values of temperature may be in doubt, the analysis should still predict 
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trends accurately. 

3.2 Definition of Temperatures 

Figure 3.1 shows a cross-section of the engine cylinder and identifies 
the various temperatures used in the analysis. Heat transfer occurs from the 
hot cylinder gases through the metal parts to the engine coolant. Three 
coolant temperatures are identified: the head coolant, the cylinder wall 
coolant and the piston coolant. The exhaust valve receives heat from the 
cylinder gases and loses it to the head by conduction through the valve seat, 
to the valve guide by conduction up the valve stem and to the exhaust port 
gases by convection. 

The temperatures labeled in Figure 3.1 are defined below. 


TCI = Piston cooling oil temperature 
TC2 = Cylinder wall coolant temperature 
TC3 = Head coolant temperature 
TWOC = Coolant side cylinder wall temperature 
TWGC = Gas side cylinder wall temperature 
TWCH = Coolant side head temperature 
TWGH = Gas side head temperature 
TWCP = Coolant side piston temperature 
TWGP = Gas side piston temperature 
TWPR = Piston rim temperature 
TWCR = Top compression ring temperature 
TWEV = Exhaust valve temperature 
TWEVG = Exhaust valve guide temperature 
TEXHP = Exhaust port gas temperature 


3.3 Network Diagram 

It is common in heat transfer problems to take advantage of the analogy 
between heat transfer through thermal resistance and direct current flow 
through electrical resistance. Steady state heat transfer relations such as 
Q = h A (T 2 - Ti) and Q = k A ^ * ^ 2 ) 
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can be written as 


where 


Q = LdL. 

y Ri 


and Q = 

*2 


R ' = r* and R * - Fa 


the temperatures Tt and T 2 represent potentials similar to voltages and the 
thermal resistances R t and R 2 are analogous to electrical resistances. One 
advantage of this approach is that it allows the heat transfer problem to be 
represented pictorial ly as a resistance network diagram. 

Figure 3.2 shows the network diagram of the heat transfer problem 
identified above. The node points in this diagram correspond to the 
temperatures identified in Figure 3.2. The resistances are defined below. 

Rj = Conductive resistance between exhaust valve and the valve guide 

R^ = Convective resistance between exhaust valve and exhaust port 
gases . 

R^ = Conductive resistance between exhaust valve and the cylinder 

head through the valve seat. 

R^ = Conductive resistance through the cylinder head. 

R = Convective resistance between head and head coolant 

D 

Rg = Conductive resistance from the gas side of piston crown to 

oil-cooled under-side of piston, 

R^ = Convective resistance on oil-cooled underside of piston 

Rg = Conductive resistance between piston and cylinder wall through the 
lower rings and the oil film. 

Rg = Conductive resistance between piston crown and piston rim. 

RlO = Conductive resistance between top compression ring and the 
cylinder wall. 

Rjl = Conductive resistance through cylinder wall. 

Rj 2 = Conductive resistance between piston rim and top compression ring. 

= Convective resistance between cylinder wall and coolant. 

^14 = ^ OIT1 ^^ ne< ^ convective and radiative resistance between cylinder 
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wall and cylinder gases. 

R = Combined convective and radiative resistance between piston crown 
Id 

and cylinder gases. 

= Combined convective and radiative resistance between head and 
lb 

cylinder gases. 

R ? = Combined convective and radiative resistance between exhaust valve 
and cylinder gases. 

These resistances can be expressed in terms of more fundamental quantities as 
follows, 

XXEVG 

K 1 " AEXVS • CEXV 
XXEVH 

K 3 _ AEXST • CEXV 

R - 1 

5 _ HCHD • AHEAD 

**7 = HOILP ♦ APIST2 
XXPCR 

K 9 " APISTR ♦ CPIST 
XXSLV 

11 _ ASLEEV • CSLEEV *12 " ARING2 • CRING 

R - 1 

k 13 " HCSL • ASLEEV 

where the heat transfer path lengths, areas, thermal conductivities, and 
convection coefficients are given below 

XXEVG = Path length from exhaust valve to valve guide 

XXEVH = Path length from exhaust valve to valve seat 

XXHED = Path length through head 

XXPIS = Path length through piston 

XXPSL = Path length from piston to cylinder wall through lower rings 
and oil film 

XXPCR = Path length from piston crown to piston rim 

XXRSL = Path length from top compression ring to cylinder wall 

XXSLV = Path length through cylinder wall 

XXPRG = Path length from piston rim to top compression ring 
AEXVS = Area of exhaust valve stem 

AEXV = Area of exhaust valve exposed to combustion chamber 

AEXST = Exhaust valve seat contact area 


R 2 “ HEXHP • AEXV 

v XXHED 

4 “ AHEAD • CHEAD 

XXPIS 

K 6 " APISH • CPIST 
XXPSL 

8 “ ASLEEV • CPIST 
XXRSL 

K 10 " ARING1 • CRING 
„ XXPRG 
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AHEAD = Cylinder head area not covered by exhaust valves 
APIST1= Cross-sectional area of cylinder 

APIST2= Surface area exposed to oil on underside of piston 
ASLEEV= Characteristic area for heat transfer through cylinder wall 
APISTR= Characteristic area between piston crown and piston rim 
ARING1= Contact area between top compression ring and cylinder wall 
ARING2= Contact area between top compression ring and piston rim 

CEXV = Thermal conductivity of exhaust valve 
CHEAD = Thermal conductivity of head 
CPIST = Thermal conductivity of piston 
CRING = Thermal conductivity of ring 
CSLEEV= Thermal conductivity of cylinder wall 

HEXHP = Convection coefficient from back side of exhaust valve to 
gases in exhaust port 

HCHD = Convection coefficient for coolant on backside of head 

HOILP = Convection coefficient for cooling oil on backside of piston 

HCSL = Convection coefficient for coolant on backside of cylinder wall 


R through R „ are combined convective and radiative resistances between 
14 17 

the time varying cylinder gas and the cylinder wall, and must be handled in a 
special way. The heat transfer to one of the combustion chamber surfaces, 
such as the head, piston, cylinder wall or exhaust valve, is given for an 
entire cycle by 

Q , = <f> h. A. (T - T ) dt 
v cyc 1 e J l l v i gas-' 

and the average continuous heat transfer rate is obtained by multiplying by 
the engine speed 

6 = »r j h i A i < T i - W dt 

The convection coefficient, h. , in this expression will be obtained from 

either the Annand or Woschni correlations already present in the Benson 

4 4 

program. The Annand correlation includes a radiation term, c (T ^ — Tg^g)’ 
this term must be included by defining a psuedo-convec t i ve coefficient as 
follows. 

« - < h c + h r> A i ' T i - T gas> 
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where h = conventional convective coefficient, 
c 

h = radiative coefficient = c [T? + T? T + T. T 2 + T 3 ] 

r L i 1 gas i gas gas J 

When using the Annand Correlation, the coefficient h^ is equal to the sum of 

h and h . 
c r 

The integral in the expression for Q can be separated and can be 
removed from its integral since it is a constant. 


Q = N T u> h, A. dt 


r i 


i i 


N r j h t A, 


T dt 
gas 


or 


Q = N © h. A. dt 
r T 1 l 


O h 

— 

1 f h i A i 


, A. T dt 
1 i gas 


dt 


By defining an effective gas temperature, T 


eg 

<t> h. A. T dt 
l i gas 


eg 


© h A A dt 


the heat transfer rate from the gas to surface i becomes 


Q = N © h. A. dt [ T. - T 1 
rjii i egj 


or 


where R = 


Q = 


T. - T 
i eg 

R 


N r j h i A i dt 

It should be noted that the effective gas temperature will be a different 
value for each surface. 


T = Effective gas temperature for exhaust valve 

eg,exh valve & 

T^g ^ ea( j = Effective gas temperature for head 

T . = Effective gas temperature for piston 

eg , p i s t 

T e g C y 2 = Effective gas temperature for cylinder wall 
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3.4 Definition of Heat Generation Terms 


In addition to the resistances and nodes shown on the network diagram in 

Figure 3.2, there are a number of heat source terms which are identified as 

either Q or Q . The Q r . terms relate to the heat generated by friction 
fric cv trie 

between the cylinder wall and the piston and rings. The Q cy terms relate to 
the heat transferred from the cylinder gases as they cool in the crevice 
volume between the piston and cylinder wall above the top ring. 


3.4.1 Friction Submodel 

The friction mean effective pressure (FMEP) due to piston/ring cylinder 

wall friction can be predicted using a correlation developed by Bishop [7]. 

4 0 2 

FMEP 6.2 X 10 r 



where p o = oil density 

U = mean pistion speed 
P 

r = compression ratio 

p 0 b 
K o P 

Re = Reynolds number = c — 

b - bore 

p o = oil viscosity 

The torque required to overcome this friction, in a two stroke engine, is 

(FMEP) V d (FMEP) b 2 S 
T f _ 2 tt “ 8 

where V, is the displacement volume and S is the stroke, 
d 

The power required to overcome the friction, and thus the rate at which energy 
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must be dissipated by heat transfer is given by 

P r = T f (j = 2 tt N T r 
f f r f 

where is the engine speed. 

The distribution of the heat transfer must be estimated either from 
measurements or more detailed analysis than that provided here. This 
information can be provided to the simulation through two defined parameters. 

CFRIC1 which is the fraction of the total piston/ring/cylinder wall 
friction which is due to the top compression ring and 

CFRIC2 which is the fraction of the heat transfer resulting from top ring 
friction which goes to the ring side. 

Using these definitions, the heat source terms for the network diagram in 
Figure 3.2 become 

Q r , = CFRIC1 • CFRIC2 * P r 

friCj f 

«fric 2 * (' - CFRIC2 > P f 
®fric - (* " C™ 01 ) ’ CFRIC2 • P f 

3.4.2 Crevice Volume Heat Transfer Submodel 

Figure 3.3 shows a cross-section of the crevice volume region of the 
piston/cylinder wall above the top ring and the dimensions which must be 
specified to characterize this volume. Gas enters and leaves this volume as 
the pressure in the cylinder increases and decreases. Heat transfer from 
these gases to the piston and top ring can significantly raise the top ring 
temperature above that due to conduction from the piston and frictional 
heating. It is possible to calculate a characteristic time for conduction as 
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t = — where L is a characteristic length and a is the thermal diffusivity of 
a 

the conducting medium. Assuming a 0.005 inch clearance between the piston and 

cylinder wall as a characteristic length and the thermal diffusivity of air at 

2000 deg F, the characteristic time becomes 

(0.005/12 ft) 2 c v ,„-5 

t = A = 6.13 X 10 sec 

2.83 X 10 ft /sec 

which is equal to 2.2 degrees of crankshaft rotation at 6000 RPM. This 
indicates that conduction will rapidly bring any air which enters the crevice 
volume into thermal equilibrium with the walls. 

The following assumptions are used in the model for crevice volume heat 
transfer . 


1. The pressure in the crevice volume is equal to the cylinder pressure. 

2. The temperature of the gas in the crevice volume is constant at a 
value determined by an area-average of the wall temperatures. 


where 


cv 


At 


P TWPR + -7— TWCR + -7— TWGC 




Af 


T = Crevice volume gas temperature 


TWPR = Piston rim temperature 
TWCR = Top compression ring temperature 
TWGC = Gas side cylinder wall temperature 
Ap = Piston area exposed to crevice volume 
Ap = Ring area exposed to crevice volume 


A^ = Cylinder area exposed to crevice volume 


At - Ap + Ar + A c 


3. Since the volume of the crevice is constant, and the temperature and 
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pressure in the crevice are known from assumptions 1 and 2, the mass of gas in 
the crevice can be calculated from the ideal gas equation 


m 


P 

cyl cv 


cv 


R T 


cv 


where P j = cylinder pressure 

V = volume of crevice 
cv 

R = Ideal gas constant for the cylinder contents 
From one time step to the next, the mass added to the crevice volume will 

be 


Am 


cv 


cv _ R T_„ ^ AP cyl^ 


CV 


where AP is the change in cylinder pressure during the time step 
cyl 

Am is the change in mass of the crevice volume 
cv 

4. Gas which enters the crevice volume is instantly cooled and this 

energy is delivered to the walls in proportion to their areas. 

(L = Am C (T - T ) 

T p v gas cv' 


and 


where 


A? 

l cv l = K Q T, 


K 

'cv n = Q T, 


} = -r 1 Qt 

cv 3 ^ 


1 * w, 2 

Q^, = energy transferred from gas to crevice volume walls. 

Am = Increment of mass added to crevice volume in one time step 


T = Cylinder gas temperature 
gas 


cv. 


cv r 


cv. 


= Fraction of heat transfer which goes to piston rim 

= Fraction of heat transfer which goes to top compression ring 

= Fraction of heat transfer which goes to cylinder wall 


Combining these equations gives the following expression for the crevice 
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volume heat transfer for one time step. 

rT 


n = f_I_ 1 v - i’ 

W tv-lj V cv T cy l \ 


AP 


cyl 


The total heat transfer for a cycle is given by summing the heat transfer at 
each time step 


Of - 


cycle 


"r-1 


V [JE5. - l| AP . 
cv I T J cyl 

L cv J 


and the average heat transfer rate is obtained by multiplying by the engine 
speed . 


\ - N r • cycle (*) V cv ~ >] iP , 


cyl 


3.5 Node Equations 

Based on the concept that at steady state the currents into any node of 

the network diagram given in Figure 3.2 must add to zero, equations can be 

written which sum the heat transfer into each node. There is one equation for 

each of the nine unknown temperatures. These equations are given below. 

TC2 - TWCC . TWGC - TWOC 


R, 


R, 


= 0 


13 11 

TCI - TWCP TWGP - TWCP 

h R 6 

TC3 - TWCH TWGH - TWCH 
Rr- R, 


= 0 


= 0 


T - TWFV 

TWEVG - TWEV TEXHP - TWEV TWGH - TWEV eg.exh valve in v 

R 1 R 2 R 3 R 17 

TWCH - TWGH TWEV - TWGH T eg.head " n 

R R + — r = 0 

K 4 K 3 K 16 

T — TWPP 

TWGC - TWGP TWPR - TWGP TWCP - TWGP eg.pist , A n 

R 8 R 9 R 6 R 15 frlc 3 = 

TWGP - TWPR TWCR - TWPR * 

R 9 R 12 CV 3 = 
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T — TWGC 

TWOC - TWGC TWCR - TWGC TWGP - TWGC A eg.cyl + + o 

R n R 10 R 8 R 14 cv 2 frlc : 

TWCR 


= 0 


TWGC - TWCR TWPR 


R 


10 


R 


12 


+ Q + Q r , =0 

cVj frlCj 


3.6 Solution Technique 

The nine node equations listed above are a system of linear algebraic 
equations and can be solved using a Gaussian elimination procedure. A 
subroutine for this purpose was obtained from Gerald [8]. The calculation is 
performed at the end of each cycle iteration providing an updated estimate of 
the nine metal temperatures. Depending on the value of the user specified 
parameter, NHEAT , four of these temperatures, the cylinder wall temperature, 
the piston crown temperature, the head temperature and the exhaust valve 
temperature, can be used in the next iteration so that the metal temperatures 
will approach their steady state values (NHEAT = 1). The other option is to 
fix these four metal temperatures at values specified by the user (NHEAT =0) . 
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4.0 INPUT DATA 


Two sets of input data must be provided to the Benson program to define 
the particular application being considered. The first data set provides 
general data describing the engine and operating conditions being 
simulated. The second data set relates specifically to the steady state heat 
transfer analysis discussed in Chapter 3.0 and contains the thermal properties 
of the engine components and the appropriate areas and path lengths. 

4. 1 Subroutine INPUT 

Below is listed the input data set required by the Benson program. 

The variable names are also provided. The current version of the program 
has these varibles entered through a BLOCK DATA statement in the main program 
and then passed through a COMMON block into the program. READ statements are 
provided in the INPUT subroutine to read an external file to obtain these 
variables but these statements have been made into comment statements. 


EQUIVD 

RPM 

PEXH 

PAIR 

TAIR 


CR 

BORE 

STROKE 

00NR0D 

FPISTA 


EVO 

EVC 

AVO 

AVC 



TWALL 

WFUEL 





ACB 

ACF 





ANNA 

FMEPM 





NPCWF 

VALEXH 

DELI 

DEL2 

DEL3 

AAN(l) 

AAN(2) 

CDE 

NTEXH 






ALPHEX(l) FEXH(l) 
ALPHEX(2) FEXH(2) 


ALPHEX(NTEXH) FEXH(NTEXH) 
VALAIR 

CDA WIDTHA 
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EQUIVD = Desired trapped equivalence ratio. If a specified fuel 
flow rate is desired, then EQUIVD sliould be set to zero 
and WFUEL set to the amount of fuel injected per cycle. 
RPM = Crankshaft speed, revolutions per minute. 

PEXH = Exhaust pressure, psia. 

PAIR = Intake manifold pressure, psia. 

TAIR = Intake manifold temperature, deg R. 

CR = Geometric Compression ratio 

BORE = Bore, inches. 

STROKE = Stroke, inches 

OONROD = Connecting rod length, inches. 

FPISTA = Ratio of piston surface area to crossectional area of 
cylinder . 

EVO = Crankangle at which exhaust valve opens, deg ATDC. 

EVC = Crankangle at which exhaust valve closes, deg ATDC 
AVO = Crankangle at which intake port opens, deg ATDC 

AVC = Crankangle at which intake port closes, deg ATDC. 

TWALL = Initial estimate of cylinder wall temperature, deg R. 
WFUEL = Initial estimate of fuel injected be cycle, unless 

EQUIVD is zero in which case WFUEL is the actual fuel 
injected per cycle, Lbm. 

ACB = Start of combustion, degrees ATDC. 

ACF = End of combustion, degrees ATDC 

ANNA = Parameter for heat transfer model 

FMEPM = Multiplier for friction model 

NPCWF = Parameter for selection of combustion model 

If NPCWF = -1 Whi tehouse-Way model used 
If NPCWF = 0 and DEL3 = 0. AVL model used 
If NPCWF = 0 and DEL3 > 0. Wiebe model used 
If NPCWF = 1 Watson model used 
If NPCWF = 2 User supplied model 
DELI = Parameter for combustion model 
DEL2 = Parameter for combustion model 
TtF.1 3 = Parameter for combustion model 

AAN(l) = Injection rate parameter for Whi tehouse-Way combustion 
model 

AAN(2) = Injection rate parameter for Whi tehouse-Way combustion 
model 

VALEXH = Parameter which determines whether an exhaust valve or 
port is used 

If VALEXH = 0. , exhaust ports are used 
If VALEXH = 10. , exhaust valves are used 
CDE = Exhaust valve/port discharge coefficient 
NTEXH = If an exhaust valve is used, NTEXH is the number of 
crankangle area combinations provided to specify the 
valve flow area 

ALPHEX(I) = Array of crank positions at which valve flow area is 
specif ied. 

FEXH(I) = Array of exhaust valve flow areas corresponding to 
ALPHEX(I) . 
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VALAIR = Parameter which determines whether an intake port or 
valve is used 

If VALAIR = 0. , intake ports are used. 

If VALAIR =10. , intake valves are used. 

CDA = Intake port/valve discharge coefficient. 

WIDTHA = If an intake port is used, this parameter specifies the 
fraction of the circumference which is taken up by port 
area. 


4.2 Subroutine INPUT2 

Below is listed the input data set required by the Steady State Heat 
Transfer Analysis. The variable names are provided below. The current 
version of the program has these variables defined within the subroutine 
INPUT2 although READ statements are provided to read these variables from an 
external file. Currently these READ statements have been made into comment 
statements . 


TWEV 

TWGH 

TWGP 

TCI 

TC2 

TC3 

NEXHV 

DEXHV 


D1 

D2 

D3 

XXEVG 

XXEVH 

XXPSL 

XXPIS 

XXPRG 

XXHED 

XXPCR 

XXRSL 

XXSLV 

AEXVS 

AEXST 


ASLEEV 

APIST2 


APISTR 

ARING1 

ARING2 

CEXV 

CHE AD 

CSLEEV 

HEXHP 

HOILP 

HCSL 


TWGC TWEVG 

D4 D5 D6 

CPIST CRING 

HCHD 


TWEV = Initial estimate of the exhaust valve temperature, deg F. 

TWGH = Initial estimate of the gas side head temperature, deg F. 

TWGP = Initial estimate of the gas side pistion temperature, deg F. 

TWGC = Initial estimate of the gas side cylinder wall 

temperature, deg F. 

TWEVG = Exhaust valve guide temperature, deg F. 

TCI = Piston cooling oil temperature, deg F. 

TC2 = Cylinder wall coolant temperature, deg F. 

TC3 = Head coolant temperature, deg F. 

NEXHV = Number of exhaust valves 


31 


DEXHV = Diameter of exhaust valves, inches 

D1 = Distance below piston surface of top ring groove, inches. 
D2 = Depth of ring groove, inches. 

D3 = Width of ring groove, inches. 

D4 = Piston/cylinder wall clearance (piston is assumed centered 
in bore), inches. 

D5 = Ring width, inches. 

D6 = Ring depth, inches. 

XXEVG = Heat transfer path length from exhaust valve to valve 
guide 

XXEVH = Heat transfer path length from exhaust valve to valve seat 

XXPSL = Heat transfer path length from piston to cylinder wall 

through lower rings and oil film 
XXPIS = Heat transfer path length through piston 
XXPRG = Heat transfer path length from piston rim to top 
compression ring 

XXHED = Heat transfer path length through head 

XXPCR = Heat transfer path length from piston crown to piston rim 

XXRSL = Heat transfer path length from top compression ring to 

cylinder wall 

XXSLV = Heat transfer path length through cylinder wall 
AEXVS = Area of exhaust valve stem 
AEXST = Exhaust valve seat contact area 

ASLEEV = Characteristic area for heat transfer through cylinder 
wall 

APIST2 = Surface area exposed to oil on underside of piston 
APISTR = Characteristic area between piston crown and piston rim 
ARING1 = Contact area between top compression ring and cylinder 
wal 1 

ARING2 = Contact area between top compression ring and piston rim 

CEXV = Thermal conductivity of exhaust valve 

CHEAD = Thermal conductivity of head 

CP I ST = Thermal conductivity of piston 

CRING = Thermal conductivity of ring 

CSLEEV = Thermal conductivity of cylinder wall 

HEXHP = Convection coefficient from back side of exhaust valve to 
gases in exhaust port 

HOILP = Convection coefficient for cooling oil on backside of 
piston 

HCSL = Convection coefficient for coolant on backside of 
cylinder wall 

HCHD = Convection coefficient for coolant on backside of head 
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5.0 OUTPUT 


The Benson program creates three output files containing the results of 
the engine simulation. The first file, written from UNIT 2. contains a 
summary of the input data and then the following results written for each 
iteration. 

1. Fuel Flow 

2. Scavenge Ratio 

3. Trapping Efficiency 

4. Scavenge Efficiency 

5. Temperature Out (Exhaust) 

6. Power 

7. True Purity 

8. Air Flow 

9. Equivalence Ratio 

10. Maximum Pressure 

11. Maximum Temperature 

12. Heat Loss (percent) 

After the program has converged a summary of the heat transfer results is 
printed showing the metal part temperatures and the fraction of the heat 
transfer going to each part of the engine. 

The second output file, UNIT 3, contains a listing of the cylinder 
pressure, temperature, specific heat ratio and volume as a function of crank 
angle for each iteration performed by the program. The last complete cycle 
of data corresponds to the cycle where the program converged. 

The final output file, UNIT 4, contains a summary of the cycle data which 
is written at the end of each iteration. It includes the trapped pressure and 
temperature, the maximum pressure and temperature, the release pressure and 
temperature, the various gas exchange efficiencies described in section 7.1 
and the indicated power. This data is printed for each iteration so the final 
listing corresponds to the listing for the converged cycle calculation. 
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6.0 SAMPLE CASES 


6. 1 Napier Nomad 

The Napier Nomad was a combined cycle diesel engine developed during the 
early 1950's as a highly efficient aircraft engine. It is the most efficient 
aircraft engine ever flown, although the fuel economy of the Nomad did not 
provide sufficient advantage over the smaller, lighter gas turbine engine. In 
this example, the Benson program was combined with a simple compressor and 
turbine model to simulate the take-off performance of the Napier Nomad engine. 
Sufficient published data of the Nomad is available to confirm the analysis. 
The Nomad was liquid cooled, and the analysis was conducted assuming both oil 
and water cooling. 

The input BLOCK DATA for this simulation is provided below. A summary of 
the results is presented in Table 6.1. A listing of the program which calls 
the Benson program including the turbomachinery is provided in Appendix 10.2. 


DATA EQUIVD.RPM.PEXH.PAIR.TAIR /0. 60, 2050. . 134.79, 149.76,896./ 
DATA CR. BORE, STROKE, C0NR0D.FPISTE.FPISTA /ll . 17,6.00,7.375, 

1 14.75,1. .1.13/ 

DATA EV0.EVC.AV0.AVC /105. ,255. , 130. ,230./ 

DATA TWALL , WFUEL /1460. .0.0001/ 

DATA ACB.ACF /335.0.10./ 

DATA ANNA , FMEPM /0. 30, 1.25/ 

DATA NPCWF.DELl ,DEL2.DEL3,AAN(1) ,AAN(2) /-I .0.01457,0.667,0. 4, 
1 0.4, 0.5/ 

DATA VALEXH , CDE , WIDTHE /0. .0.85,0.20/ 

DATA VALAIR , CDA . WIDTHA /0. .0.85.0.36/ 
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Table 6.1 


Napier Nomad 


Specifications : 

Two Stroke, Loop Scavenged 

12 Cylinders Bore: 6.00 inches Stroke: 7.375 inches 

Effective Compression Ratio: 8.00 

At Takeoff Power: Speed: 2050 rpm Mean Piston Speed: 2520 ft/min 

Inlet Manifold Pressure: 89 psia 
Inlet Manifold Temperature: 475 F 
Compressor Efficiency: 0.85 
Turbine Efficiency: 0.84 
Under Takeoff Conditions 

Published Data Benson Program Benson Program 

(Water Cooled) (Oil Cooled) 


BMEP 

Peak Pressure 
Air Flow 


Compressor Power 
Total Shaf t Power 
Overall BSFC 
Percent Heat Loss 
Head Temp 
Piston Temp 
Sleeve Temp 


205 psi 
2200 psi 
13 lb/sec 
1222 F 


2250 hp 
1840 hp 
3070 hp 
0.345 lb/hp-hr 
NA 
NA 
NA 
NA 


205. 1 psi 
2186 psi 
13.03 lb/sec 
1162 F 
2654 hp 
2374 hp 
1858 hp 
3169 hp 
0.339 lb/hp-hr 
13.8 % 

782 F 
834 F 
791 F 


205.4 psi 
2205 psi 
13.00 lb/sec 
1199 F 
2658 hp 
2424 hp 
1854 hp 
3228 hp 

0.329 lb/hp-hr 

10.2 % 

1498 F 
1008 F 
1224 F 


Turbine Inlet Temp 
Diesel Power (12 cyl) 2660 hp 
Turbine Power 
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6.2 Garrett OCE 


The Garrett Turbine Engine Company, under Army contract, used the Benson 
program to simulate a compound cycle engine suitable for use in a light 
helicopter. The engine proposed by Garrett is a 6 cylinder, two stroke diesel 
with two expansion turbines on the diesel exhaust. The first turbine drives 
the compressor for the diesel intake air and the second is a power turbine 
which is connected directly to the engine’s output shaft. As in the above 
example, the Benson program was combined with a simple compressor and turbine 
model to simulate the Garrett design for a Compound Cycle Engine. Data for 
the predicted performance was supplied by Garrett in the final contract 
report. This data is compared below in Table 6.2 with the results of the 
simulation for both water and oil cooling. The input data set is also 
provided below. A listing of the main program used to simulate the Garrett 
engine is provided in Appendix 10.3. 


DATA EQUIVD , RPM , PEXH , PAIR , TAIR /0. 68. 6122. . 134.79, 149.76,896./ 

DATA CR, BORE, STROKE. OONROD.FPISTE.FPISTA /9. 1712,3. 101 ,2.940, 

1 7.168,1. ,1.13/ 

DATA EVO , EVC , AVO , AVC /90. .239. ,126. ,234./ 

DATA TWALL , WFUEL /I 460. ,0.000145332/ 

DATA ACB.ACF /345..10.5/ 

DATA ANNA , FMEPM /0. 36. 1.0/ 

DATA NPCWF.DELl , DEL2 , DEL3 , AAN ( 1 ) ,AAN(2) /-I ,0.01457,0.667,0.4, 

1 0. 4,0.5/ 

DATA VALEXH , CDE , NTEXH /10. ,0.8864, 13/ 

DATA VALAIR , CDA , WIDTHA /0. ,0.8,0.5542/ 

DATA ALPHEX/0. 00. 10. 97. 25. 00, 31. 61, 41. 00, 52. 25. 72. 89, 90. 95. 104. 00, 
1 111.59,119.00, 132.23, 149.00,37*0./ 

DATA FEXH/0 .0,0. 00093264 , 0 . 00364 , 0 . 00591222 ,0.00899 ,0.0106201 , 

1 0.01175263,0.01111131,0.00923,0.00717895.0.00461,0.00181743.0. , 

2 37*0./ 
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Table 6.2 


Garrett QCE 

Specifications : 

Two Stroke, Uni flow Scavenged 

6 Cylinders Bore: 3.101 inches Stroke: 2.940 inches 

Effective Compression Ratio: 7.50 

At Rated Power: Inlet Manifold Pressure: 149.79 psia 

Inlet Manifold Temperature: 436 F 
Compressor Efficiency: 0.79 
Turbine Efficiency: 0.845, 0.848 
Rated Power Conditions 


Published Data 

Benson Program 
(Water Cooled) 

Benson Program 
(Oil Cooled) 

BMEP 

400 psi 

386.8 psi 

386.5 psi 

Peak Pressure 

3364 psi 

3401 .1 psi 

3429.2 psi 

Air Flow 

2.349 lb/sec 

2.316 lb/sec 

2.301 lb/sec 

Turbine Inlet Temp 

1752 F 

1738 F 

1785 F 

Diesel Power (6 cyl) 824 hp 

796.0 hp 

795.5 hp 

Turbine Power 

226 hp 

233.1 hp 

248.6 hp 

Fan Power 

51 hp 

51 hp 

51 hp 

Total Shaft Power 

999 hp 

978.1 hp 

993.1 hp 

Overall BSFC 0. 

328 lb/hp-hr 

0.332 lb/hp-hr 

0 . 323 lb/hp-hr 

Percent Heat Loss 

8.46 % 

8.59 % 

5.74 % 

Head Temp 

1244.3 F 

1299.6 F 

2078.3 F 

Piston Temp 

1279.3 F 

1375.8 F 

1571.3 F 

Sleeve Temp 

418.3 F 

1344.8 F 

1792.2 F 
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7.0 DOCUMENT AT I OF SUBROUTINES 


7. 1 Subroutine BENSON 

A simulation using the Benson program is initiated by a call to 
subroutine BENSON. There is no argument list for this subroutine. All data 
is passed into the subroutine through input files or COMMON blocks and all 
output data is written to the three output files discussed in Chapter 5.0. 

Subroutine BENSON calls the input subroutines, the initialization 
subroutine, and then organizes the calls to CYLNDR and POWER to integrate the 
energy balance throughout the cycle. Examples of how BENSON is called from a 
main program are given in the Appendix for the sample cases of Chapter 6.0. 

7.2 Subroutine BAL 

This subroutine can be used for two purposes: 

1. To find the enthalpy of combustion products given the fuel-air 
ratio and the temperature. 

2. To find a final temperature of combustion products given an 
initial temperature and a rise in the enthalpy of the products. 

For the combustion equation 

A • CH n + [0.21 0 2 + 0.79 N 2 ] -♦ Xj N g + 0 2 + ^ C0 2 + X 4 H 2 0 

the number of moles of products, per mole of air. are given by the following 
equations . 

X x = 0.79 

^ = 0.21 - Xj - i X 4 
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lb fuel # lb C # lb air # mole C _ mole C 

*3 - lb air lb fuel mole air lb C " mole air 

= FZA • CARBON • 28.85 / 12.01 

_ lb fuel t lb H 2 # lb air # mole H 2 _ mole H 2 

*4 - lb air lb fuel mole air lb H 2 - mole air 

= FZA • (1 - CARBON) • 28.85 / 2.016 

These quantities can be used to calculate mole fractions by summing the 
total number of moles and dividing each term by the total. 

To accomplish the first purpose of this subroutine the user specifies the 
temperature and fuel/air ratio of the mixture as well as the mass fraction of 
carbon in the fuel and the subroutine calculates the mole fractions of the 
product gas composition and then evaluates fourth order polynomial curve fits 
to JANAF table data to obtain the enthalpy of the gas mixture at the specified 
temperature [1]. 

For the case where the subroutine is intended to locate a temperature 
given an initial temperature and an increment of enthalpy, the subroutine 
first evaluates the enthalpy at the specified initial temperature and adds the 
enthalpy increment to obtain the enthalpy of the unknown state. Then the 
subroutine guesses a temperature and tests whether the enthalpy at this 
temperature is the same as the unknown. If not, the program adjusts the 
temperature and checks again. Newton’s method is used to speed convergence. 
When the enthalpy matches within 0.01%, the subroutine returns to the calling 
program. 

The argument list for subroutine BAL is 

SUBROUTINE BAL(CARB0N.FZA,T1 ,DELH,T2,GAM,H3) 
where CARBON is the mass fraction of carbon in the fuel 
FZA is the fuel/air mass ratio 
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T1 is the initial temperature when the subroutine is being used to 
locate a second unknown temperature. 

DELH is the enthalpy increment. 

T2 is the temperature at which the enthalpy is desired or is the 
unknown temperature. 

GAM is the ratio of specific heats evaluated at T2. 

H3 is the enthalpy at T2. 

7.3 Subroutine CYLNDR 

This subroutine integrates the energy equation described earlier for the 
gas exchange portion of the cycle. The subroutine checks to see if the valves 
or ports are open and, if so, integrates the energy equation. If the valves 
are closed, then CYLNDR calls POWER to integrate the energy equation for the 
compression-expansion portion of the cycle. When the subroutine is called the 
first time step after the gas exchange has been completed, the subroutine 
totals the mass added and removed from the cylinder and calculates the various 
gas exchange efficiencies. These efficiencies are defined below. 


^ . mass of delivered air retained in cylinder 

Scavenging Efficiency. mass of trapped cylinder cteriT 


Purity: 


mass of air in trapped cylinder charge 
total mass of trapped cylinder charge 


Charging Efficiency based on Trapped Condition: 

mass of delivered air retained in cylinder 
Trapped volume x trapped density 
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Charging Efficiency based on Inlet conditions and trapped volume: 

mass of delivered air retained in cylinder 
Trapped volume x inlet density 


Charging Efficiency based on inlet conditions and swept volume: 

mass of delivered air retained in cylinder 
Swept volume x Inlet density 


Trapping Efficiency: 


mass of delivered air retained in cylinder 
total mass of air delivered to cylinder 


Scavenge Ratio based on swept volume: 

total mass of air delivered to cylinder 
Swept volume x inlet density 


Scavenge Ratio based on trapped volume: 

total mass of air delivered to cylinder 
Trapped volume x inlet density 


Relative Charging Efficiency based on trapped volume: 

total mass of trapped cylinder charge 
Trapped volume x inlet density 


Relative Charging Efficiency based on swept volume: 

total mass of trapped cylinder charge 
Swept volume x inlet density 


7.4 Subroutine CYLVOL 


This subroutine takes the piston position and the piston velocity 
provided by STROK and calculates the cylinder volume and the rate of change of 
volume. 
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The volume at some crank position is equal to the sum of the volume 
displaced by the piston plus the displacement volume. The volume displaced by 
the piston is equal to the piston area times the distance between the top of 
the piston and its position at top dead center which is XSTROK as calculated 
by the subroutine STROK. The clearance volume can be expressed in terms of 



AREA is the piston area 

r is the compression ratio 

Thus, the volume in the cylinder can be expresed as 

q . apfA 

V = AREA • XSTROK + 2 

or 

V = AREA • [XSTROK + - - j -y ] 

The rate of change of the volume is equal to the derivative of this expression 
or in terms of DSTROK provided by the subroutine STROK, 

^ = AREA • DSTROK 

uu 

The argument list follows. 

SUBROUTINE CYLVOL ( FCYLG , XSTROK , STROKE . CRR , DSTROK , VCYLG , DVCYLG ) 
where FCYLG is the cross-sectional area of the cylinder 
XSTROK is the piston position 
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STROKE is the engine stroke 
CRR is the compression ratio 
DSTROK is the piston velocity 
VCTLG is the volume 

DVCLYG is the rate of change of the volume 

7.5 Subroutine ELIH 

This subroutine is used to solve the system of nine linear algebraic 
equations developed by subroutine HTRANS for the metal part temperatures in 
the engine. The subroutine was copied from Gerald [8] and uses Gaussian 
elimination with partial pivoting. The argument list is as follows. 

SUBROUTINE ELIM(AB.N,NP.NDIM) 

where AB = Coefficient matrix augmented with vectors consisting of the 

right hand sides of the system equations (may be more than one 
set) . 

N = Number of equations. 

NP = Total number of columns in the augmented matrix AB. 

NDIM = First dimension of matrix AB in the calling program. 

The solution vector(s) are returned in the augmentation columns of AB. 

7.6 Subroutine FCLOSE 

This subroutine closes the input/output files and states that the files 
created should be kept after completion of the program. 
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7.7 Subroutine FOPEN 


This subroutine opens the input and output files which the Benson program 
reads from and writes to. Unit 1 is the input file and it is named 
DIESEL1.INP and units 2 through 5 are output files with names that can be 
specified by the user in the calling program, the user defined names are 
passed to FOPEN as 16 character names in the COMMON block named OPEN. 

7.8 Subroutine HEADER 

This subroutine writes titles and column headings to the output files and 
also writes out selected values of the input parameters so the particular run 
can be identified. The subroutine has no argument list and all data is passed 
into the subroutine through COMMON blocks /GEN/ and /CYL/. 

7.9 Subroutine HTRANS 

This subroutine executes the procedure described in Section 2.2 for 
calculating the engine part temperatures. The subroutine calculates the 
resistances described in that section and assembles the coefficient matrix for 
the system of algebraic equations for solution by subroutine ELIM. HTRANS 
also calculates the friction mean effective pressure due to the piston/ring 
friction so that the heat transfer terms associated with the dissipation of 
this friction can be added to the coefficient matrix. The terms for the 
crevice volume heat transfer are evaluated in subroutine POWER but are added 
to the coefficient matrix in HTRANS. 

After solving the system equations by calling subroutine ELIM, the 
various heat transfer flow paths in the network diagram are summed to provide 
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a summary of the distribution of the heat transfer throughout the cylinder 
area. This data is provided to the main program through OOMMON/HEAT1/ and 
C0MM0N/HEAT2/ and may be printed if desired. 

7. 10 Subroutine INITLZ 

This subroutine is called early in the execution of the Benson program to 
perform a variety of tasks necessary to starting the calculations. If 
necessary, input variables are converted to the proper dimensions and then 
many of the variables are nondimensionalized. Variables used in summing 
calculations are assigned initial values of zero. The values of the 
polynomial coefficients used to calculate gas properties are assigned in this 
subroutine based on curve fits to the JANAF table data. The subroutine has no 
argument list and all data is passed into the subroutine through COMMON blocks 
/XXX/, /GEN/ and /CYL/. 

7.11 Subroutine INPUT 

This subroutine is used to supply the Benson program with input data 
about the specific engine application being simulated. The subroutine 
contains read statements for the essential data, however, these statements 
have been converted to comment statements and the data is either present as 
assignment statements in the subroutine or is assigned in the main calling 
program and then passed into the subroutine through 00MM0N/BNSN/ . This was 
done so that the Benson program could be repeatedly called as a subroutine 
itself without having to reread the input data file for every run. INPUT has 
no argument list and all variables are passed through to the other subroutines 
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using COMMON blocks. 


7.12 Subroutine INPUT2 

This subroutine is used to supply the Benson program with input data for 
the detailed heat transfer calculations. The subroutine contains READ 
statements for the data, however, these statements have been converted to 
comment statements and the data is assigned directly to the variables by 
executable statements within the subroutine. To change the values of these 
parameters the subroutine must be edited and then recompiled. The format of 
the input data file and the variable name list are given in Section 3.2. 


7. 13 Subroutine MASSFL 

Flow through the intake and exhaust ports and valves can be modeled as an 
adiabatic quasi-steady flow process. Potential energy effects are generally 
assumed negligible. With these assumptions, the energy balance for a control 
volume takes the following form. 


y2 y2 

h l + 2 = h 2 + *2 

where state 1 designates the upstream condition and state 2 is the condition 
at the minimum area section of the flow passage. Since state 1 is normally 
taken to be a stagnation condition, V! can be assumed to be zero. Thus, 

V 2 = [2(h 1 - h 2 )] 1/2 


or for constant specific heat 

V 2 = [2 C p (T, - T 2 )] 1/2 

Since the process is assumed isentropic, again for constant specific heats. 
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Tf-1 


II 

t 2 


Pi 

[P2. 


where t is the ratio of specific heats, C^/C^. 

Substituting this expression into the equation for V 2 gives 


V 2 = 


2 C p T 2 


[#-■] 


1 1/2 


if this velocity is divided by the speed of sound at state 2, a 2 = viRT , 
then a Mach number is obtained 


U = — ^ = 
m a 2 


TT-1 


1 1/2 


W [ft] ’ - 


if this Mach number is greater than or equal to l. f then the flow through the 
valve or port is choked and the velocity at section 2 will be sonic. If the 
Mach number is less than 1. then the flow is not choked and the velocity at 
section 2 will be given by the above expression. 

The mass flow rate through the valve or port can be obtained from the 
continuity equation 


m = p A V 

When this equation is applied at section 2 the following expression for mass 
flow through the valve or port is obtained. 


m — P 2 A 2 3-2 


TT-1 


»r->] 


1/2 


In the Benson program, the expression for mass flow rate is 

non-dimensionalized using several defined reference quaint i ties. A list of 

these quantities follows. 

T = Reference Temperature 

P r = Reference Pressure 
ref 
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= Reference Volume 


ref 

F r = Reference Area 
ref 

X . = Reference length 
ref 

D r = Reference diameter 
ref 

a re j. = Reference speed of sound, a re f = V'tRT 
= Reference mass, m 


P .V f tP .V . 
ref ref ref ref 


m 


ref 


ref 


R T 


ref 


ref 


Using these reference quantities, which are assigned internally in the 
program, the mass flow equation can be rearranged to obtain the form used in 


the Benson program. 


■ * *2 A 2 a 2 U m 


but 


P 2 = Pj 


P 

PJ 


and 


a 2 = V,rRT 2 


m = p 


1 


P 

PJ 


A 2'^ 5T I U m 


Now, from the equation defining the Mach number, 


p*_ 

Pi " 


1 

t m 


TT 

1-T 


So the mass flow can be arranged as 


** P > A 2 U m 


m = 


1 + 


rlu 2 1 

t m J 


Tr+1 

2(-r-l) 


The non-dimensional mass flow rate, WFLOW, is equal to 

- rJ ™ 

— rof J L rpf J 


m 

WFLOW = u = L ‘ref J L "ref J 

m r V r l 
ref ref 


P ref a ref J 


AO t 


ref 


where AO = — 1 — 
a ref 
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t U 


and 


FM = 


1 + 


TT-1 


m 


U 2 


nr+1 

2(t-1) 


2 m 

if the flow is choked, the equation is the same but FM becomes 

nr+1 

rr o ' 

FM = -T 


1/2 


[far 1 : 

The argument for subroutine MASSFL follows. 


SUBROUTINE MASSFL (FHOLE , RUP , RDP , AO , GR , G , WFLOW) 


where FHOLE is the non-dimensional flow area, 


ref 


RUP is the upstream pressure divided by the reference pressure 
RDP is the downstream pressure divided by the reference pressure 
AO is the ratio of the speed of sound at the upstream temperature 
to the reference speed of sound. 

GR is the reference value of t, the ratio of specific heats 

G is t, the ratio of specific heats 

WFLOW is the non-dimensional mass flow rate 


7.14 Subroutine POWER 

This subroutine integrates the energy balance equation for the 
compression-expansion portion of the cycle. The techniques used are the same 
as those described in detail in section 2.1 of this report. 

Power is called by CYLNDR after the valves close indicating the start of 
the compression process. POWER calculates the gas properties at the start and 
end of each time step to obtain the internal energy and calculates the work 
and heat transfer. POWER also contains the combustion calculations resulting 
from the four combustion models. The amount of fuel burned is calculated and 
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multiplied by the constant volume heating value to obtain an effective heat 
release. This heat release is included in the energy balance. At the end of 
the expansion process, the indicated power and efficiency is calculated and 
printed to a file. 

7. 15 Subroutine SPLFIT 

This subroutine is used to interpolate between discrete data provided to 
define a valve flow area curve. The calling program provides the X and Y 
arrays of discrete data and a value of X at which an interpolated Y is 
required and the subroutine fits a spline curve to the data and evaluates the 
spline at X. The relevant argument list follows. 

SUBROUTINE SPLFIT (0,0,NP,0,0. .XKX.O.FY.O.X.Y, ITRP) 
where NP is the number of points describing the curve 

XKX(I) is the NP X-values defining the curve (25 points max). These 
points must be in ascending order 

FY(I) is the array of NP Y-values defining the curve (corresponding 
to the X-values) 

X is the X-value where the function value is required 
Y is the calculate function value at the desired X-value 
ITRP is an error indicator. 

7. 16 Subroutine SREA 

This subroutine is used to calculate the surface area of the portion of 
the combustion chamber exposed to the combustion gases at a given crank 
position. This area will be the sum of the areas of the head surface, the 
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piston surface and the instantaneous exposed cylinder wall surface. Each of 
these areas is calculated as follows. 

Head Area = FPISTE | (BORE) 2 
Piston Area = FPISTA | (BORE) 2 
Cylinder Wall Area = tt (BORE) L 

where FPISTE = the ratio of the head surface area to the cylinder 
cross-section 

FPISTA = the ratio of the piston surface area to the cylinder 
cross-section 

BORE = the cylinder diameter 

L = the distance between the current position of the piston and its 
position at top dead center 

7 r 2 

Since the volume at any time is approximately equal to ^ (BORE) L 
it is possible to express the cylinder wall surface area in terms of the 
volume as 

4 4V 

Cylinder wall area = tt (BORE) = 

t r (BORE) BORE 

Using this relation, the total surface area can be expressed as 

AREA = + (FPISTE + FPISTA) £ (BORE) 2 

dUKE 4 

The subroutine SREA evaluates this equation for a given set of input 
arguments. The argument list follows. 

SUBROUTINE SREA (SAREA,FCYLL,DCTLL,XR, DR, FPISTA, FPISTE) 
where SAREA is the total surface area 
FCYLL is equal to J (BORE) 2 
DCYLL is the Bore 

XR and DR are reference dimensions (assumed = 1 . in the BENSON program) 
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7.17 Subroutine STROK 


This subroutine calculates the piston position and the rate of change of 
piston position (piston velocity) for a specified engine geometry and at a 
specified crank angle. Ferguson [9] gives the following expression for the 
instantaneous cylinder volume for a reciprocating engine. 


V = V 


cl 


1 + 


|l - cos 0 + i{l-(l-fc 2 


where V , is the clearance volume 
cl 

r is the compression ratio 

2 L 

e is a geometric ratio = — g— 

L is the connecting rod length 

S is the stroke 

0 is the crank angle 

By defining the quantity FNN 



FNN 


sinV 


1/2 


J 


= tl 2 ' 

and recognizing that the compression ratio is equal to y 


V . + V . 
d cl 


so that 


cl 


V , 

( r -l) _ — it can be shown that the volume can be expressed as 


cl 


V = V 


•1 + 5T f 1 " cos 6 + e ~ ¥m \ 


or 

2 ^ 
V = V cl + T/4 2 b s [l - cos 0 + \ - FNNj 

where b = bore 

V, = displacement volume 
d 

The subroutine STROK calculates the quantity 
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X = | 1 - cos 0 + i - FNnJ 

which is the distance between the top of the piston at the specified crank 
angle and its position at top dead center. Thus, at TDC, when 0=0, the 
value of X = 0. This value is provided in non-dimensional form by dividing by 
a reference length, XREF 

xs™* - m? * HM 1 - cos e * i - ™] 


The relation for the piston position can be differentiated with time to 
provide an expression for the piston velocity. 


d[X] S fdfcos 0) _ d(FNN)l 

dt " 2 [ dt dt J 



sin 0 + 


sin 20l 

2-FNN J 


d0 

dt 



sin 0 + 


sin 20 
2 FNN . 


• 2v • REVENG 


where REVENG is the engine speed in rev/sec. 

The piston speed can be non-dimensional i zed by introducing a non-dimensional 
time parameter, ZREF = where REVREF is a reference engine 

speed and AREF is a reference speed of sound. 

k™ 0 * - AHEFpdT 1 ] = m ZREF [l$ii] [sab] [ sln 6 * l^sr] 

The non-dimensional piston position and velocity supplied by subroutine STROK 
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are used by other routines in the Benson program to calculate such quantities 
as the cylinder volume and cylinder surface area. 

The argument list follows. 

SUBROUT I NE STROK ( XSTROK . DSTROK , THETAE . STROKE . CONROD . RPS , ZR ) 
where XSTROK is the nondimensional piston position returned by STROK 
DSTROK is the nondimensional piston velocity returned by STROK 
THETAE is the crank position 
STROKE is the stroke normalized by XREF 
OONROD is the connecting rod length normalized by XREF 
RPS is the rotational speed of the engine, normalized by 
dividing by the reference engine speed 
ZR is the nondimensional time parameter (ZREF) 


7 . 18 Subroutine UNITP 

This subroutine is used to change the units of a pressure supplied as an 
input. The argument list follows. 

SUBROUTINE UNITP (IA, IB.PP.PPB) 

where IA and IB are used to designate the type of unit conversion desired, PP 
is the pressure which will be returned after having been multiplied by the 
appropriate conversion factor and PPB is the barometric pressure. PPB is 
required for conversions from gage pressure to absolute pressure. The values 
of IA and IB should be selected from the following choices. 


IA IB 

1 1 

1 2 

2 1 

2 2 

* 1.2 1 

* 1,2 2 


Conversion 

do nothing 
do nothing 

psia -» bars (absolute) 
bars (absolute) -» psia 
psig -* bars (absolute) 
bars (gage) -» psia 
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7. 19 Subroutine UN ITT 


This subroutine is used to change the units of a temperature supplied as 
an input. The argument list follows. 

SUBROUTINE UNITT (IA.IB.TT) 


where IA and IB are used to designate the type of unit conversion desired and 
TT is the temperature which will be returned after having been multiplied by 
the appropriate conversion factor. The values of IA and IB should be selected 
from the following choices. 


IA 

IB 

Conversion 

1 

2 

None 

2 

2 

oK -» oC 

3 

2 

oK -> oR 

/ 1,2,3 

2 

oK -» oF 

1 

* 2 

None 

2 

* 2 

oC oR 

3 

* 2 

°R -* oK 

* 1.2,3 

* 2 

oF -* °K 


7.20 Subroutine UNITW 

This subroutine is used to change the units of a mass supplied as an 
argument. The argument list follows. 

SUBROUTINE UNITW (IA.IB.WW) 

where IA and IB are used to designate the type of unit conversion desired and 
is the mass which will be returned after having been multiplied by the 
appropriate conversion factor. The values of IA and IB should be selected 
from the following choices. 
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IA 


IB 


Conversion 


* 

* 


1 

1 

1 

1 


2 


* 2 
2 


* 2 


None 

None 

kg -* lbm 
lbm -* kg 
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8.0 CONCLUSION 


8. 1 Suggestions for Future Work 

The major weakness of thermodynamic modeling diesel engines is the lack 
of a reliable and simple way to predict the heat transfer and heat release 
rates. Correlations are available for both these quantities and they are 
available in the Benson program, however, these correlations were developed 
under conditions very different from those present in the proposed Compound 
Cycle Engine. Experimental measurements at speeds above 6000 RPM and boost 
pressures above 100 psia are not currently available in the technical 
literature. These measurements must be made and used to verify the 
correlations before they can be used reliably at these conditions. 

The current combustion models in the Benson program are of two types. 

The Watson model , the AVL model and the Wiebe model are simple functional 
models and cannot be used to understand the details of the diesel combustion. 
The Whi tehouse-Way model is an attempt to include some of the physics of the 
fuel injection and combustion processes in a model and so it has the potential 
to provide the predictive information required for design. However, the 
Whi tehouse-Way model is based on an out-of-date understanding of diesel 
combustion. This model treats the spray from a droplet perspective whereas 
more modern approaches recognize that the droplets vaporize very quickly after 
entering the combustion chamber (especially in a highly rated engine) and the 
subsequent combustion has the characteristics of a gaseous turbulent diffusion 
flame. The more successful modern diesel combustion models are based on this 
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approach. A more modern combustion model should be added to the Benson 
program allowing the effects of parameters such as fuel injection pressure, 
injection duration, nozzle diameter, air swirl to be investigated. 

Another useful feature that should be included in a model of the complete 
Compound Cycle Engine is second law thermodynamic analysis. Second law 
analysis is based on tracking the utilization and destruction of availability 
through the engine. Availability is defined as ’’the maximum amount of work 
that could be obtained by bringing a system into equilibrium with its 
surroundings.’’ Availability enters the engine with the fuel and either leaves 
as work or is progressively destroyed by irreversible processes inside the 
engine. In contrast to conventional thermodynamic analysis which tracks 
energy through the system and cannot tell whether the energy could be 
profitably converted to work, second law analysis allows the sources of 
irreversibility in the engine to be identified and quantified so that their 
effect on the fuel efficiency of the engine can be determined. 

8.2 Summary 

This report documents the Benson Diesel Engine Simulation Program. This 
program is currently configured to simulate a loop or uniflow scavenged two 
stroke engine. The program can provide information about the effects of 
combustion timing, valve timing, heat release rate and other engine design and 
operating parameters. The program has been extended to include a more 
detailed heat transfer model that allows the prediction of metal part 
temperatures. 
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10.0 APPENDIX 


10.1 Variable Name List 


ACB = Start of Combustion 

ACF = End of Combustion 

AIR = An indicator of whether the intake valve/port is open or 
closed during the cycle (=0 when closed, =10 when open) 
AIRVMX = Maximum flow area of intake valves/ports 
ALPAIR = Array of crank positions corresponding to the intake 
valve flow areas given in FAIR, degrees. 

ALPHAT = Cummulative total crankangle during integration. 

ALPHEX = Array of crank positions corresponding to the exhaust 
valve flow areas given in FEXH, degrees. 

ANGEND = Maximum number of crank degrees allowed before 
integration is terminated. 

AAN(l) = Fuel injection parameters for Whi tehouse-Way model. 

AAN(2) = Fuel injection parameters for Whi tehouse-Way model. 

AVC = Intake valve/port closing angle, degrees. 

AV0 = Intake valve/port opening angle, degrees. 

AREF = Reference value for speed of sound (based on TREF). 

ANNA, ANNB, ANNC = Parameter for Annand or Woschni Heat Transfer 
mode 1 . 

BETA = Parameter in Watson combustion model. 

CARBON = Mass fraction of carbon in fuel. 

CALVAL = Constant volume heat of combustion with water present as 
vapor in the products. 

CD1, CD2 = Parameters in Watson combustion model. 

CDA = Discharge Coefficient for intake valve/port. 

CDE = Discharge Coefficient for exhaust valve/port. 

CEFFYI = Charging efficiency based on inlet conditions 
and trapped volume. 

CHEFFY = Charging efficiency based on trapped conditions. 

00NDY = Thermal Conductivity of cylinder contents. 

CYCLE = 2 for a two stroke engine, = 4 for a four stroke engine. 
OOEFFA, C0EFFB, 00EFFC, C0EFFD = Coefficient for polynomial fit 
to JANAF tables. 

CPI, CP2 = Parameter in Watson Combustion model. 

CPT = Specific heat at constant pressure. 

CVT = Specific heat at constant volume. 

DELI, D EI ,2 . DEL3 = Parameter in either Wiebe or Whi tehouse-Way 
combustion models. 

DHEAT = Heat transfer rate from cylinder during gas exchange. 
DREF = Reference diameter. 


60 



DZ = Nondimensional time increment. 

DQF = Heat addition to cylinder due to combustion. 

DQCR = Heat transfer from cylinder due to radiation. 

DQCY = Heat transfer from cylinder due to convection. 

ETI = Internal energy of gases in cylinder at start of 
time step. 

ETF = Internal energy of gases in cylinder at end of time 
step. 

EVC = Exhaust valve/port opening angle, degrees. 

EVO = Exhaust valve/port closing angle, degrees. 

EXH = An indicator of whether the exhaust valve/port is open 
or closed during the cycle. (=0 when closed, =10 when 
open) . 

EXHVMX = Maximum flow area of exhaust valves/ports. 

FCYL = Crossectional area of cylinder. 

FMEP1, FMEP2, FMEP3 = Parameters in FMEP correlation. 

FMEPM = Multiplier for FMEP correlation. 

FPISTA = Ratio of surface area of piston top to crossectional 
area of cylinder. 

FPISTE = Ratio of surface area of head exposed to combustion 
chamber to crossectional area of cylinder. 

FAIR = Array specifying the geometric intake valve flow area. 

FEXH = Array specifying the geometric exhaust valve flow area. 

FREF = A reference area used to normalize areas in the program. 
FZA = Fuel/air mass ratio. 

GA = Ratio of specific heats for air. 

GE = Ratio of Specific heats for exhaust products. 

GREF = Reference ratio of specific heats. 

HEATG = Total heat transfer during gas exchange. 

HEATRF = Heat released due to combustion 

HC1 ... HC8 = Integral of the convective coefficient times the 
area for steady state heat transfer analysis. 

HCT1 ... HCT8 = Integral of the convective coefficient times 
the gas temperature times the area for steady 
state heat transfer analysis. 

HTA = Enthalpy of air entering cylinder 

HTF = Enthalpy of cylinder contents at end of time step. 

NPCWF = Parameter which determines which combustion model is 
used in the program 

when NPCWF. LT.O Whi tehouse-Way Model 

when NPCWF. EQ.O and DEL3.EQ.0 AVL Model 
when NPCWF. EQ.O and DEL3.GT.0 Wiebe Model 
when NPCWF. CT.O Watson Model 

NTAIR = Number of points provided in ALPAIR and FAIR to specify 
the intake valve flow area. FAIR(NTAIR), ALPAIR(NTAIR) 
NTEXH = Number of points provided in ALPHEX and FEXH to specify 
the exhaust valve flow area. FEXH ( NTEXH ) , ALPHEX (NTEXH) 
PBARAB = Barometric pressure. 

PI =3.141592654 

P02 = Partial pressure of oxygen. 

PREF = Reference pressure. 



PORTED = 1 for loop scavenged, = 2 for uni flow scavenged. 

PVEL = Mean piston speed. 

PURITY = Purity. Ratio of mass of air in trapped cylinder charge 
to total mass of trapped cylinder charge. 

RELCHG = Relative charging efficiency based on trapped volume. 
RELCHS = Relative charging efficiency based on swept volume. 
REVENG = Engine speed, rev/sec. 

REVREF = Reference engine speed, rev/sec. 

REY = Reynolds number of cylinder contents. 

SCAEFF = Scavenging Efficiency 

9CASW = Scavenge ratio based on swept volume. 

9CAVR = Scavenge ratio based on trapped volume. 

STOIC = Stoichiometric fuel/air ratio. 

SURFC = Surface area of cylinder wall exposed to combustion. 
SURFH = Surface area of head exposed to combustion chamber. 

SURFP = Surface area of piston. 

TAIR = Temperature of air entering cylinder through intake. 

TEXH = Exhaust temperature (initial estimate specified as 
input) 

TREF = Reference temperature. 

TREFF = Trapping Efficiency. 

TWALL = Cylinder wall temperature. 

VALAIR = 0 if engine has intake ports, = 10 if engine has intake 
valves . 

VALEXH = 0 if engine has exhaust ports, = 10 if engine has 
exhaust valves. 

VCYL = Cylinder volume. 

VEFFYI = Charging efficiency based on inlet conditions and 
swept volume. 

VISCTY = Viscosity of cylinder contents. 

WABU = Mass of air required to burn fuel under stoichiometric 
conditions . 

WAR = Water to air ratio. 

WCA = Air retained in the cylinder. 

WCAIR = Total mass of air trapped in cylinder. 

WCHG = Mass of air plus residual in cylinder per cycle. 

WOO = Total mass flow rate through exhaust valve. 

WCYL = Total mass in cylinder. 

WCYLT = Trapped mass in cylinder. 

WIDTHA = Fraction of cylinder circumference allocated for intake 
ports . 

WIDTHE = Fraction of cylinder circumference allocated for exhaust 
ports . 

WFUEL = Mass of fuel injected per cycle. 

WMW = Molecular weight. 

WN(I) = Number of moles of species I. 

WNA = Number of moles of air trapped in cylinder. 

WNI(I) = Number of moles of each species present at start of time 

step. 

WNF(I) = Number of moles of each species present at end of time 
step. 
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WNP = Increase in number of moles in cylinder due to 
combustion (moles air - moles products). 

WNR = Number of moles of residual trapped in cylinder. 

WNTF = Total number of moles of gases in cylinder. 

WRES = Mass of residual gases in cylinder. 

WTOT = Total mass in cylinder, air+residual+fuel 



10.2 Napier Nomad program listing 


C 

C 

c 


c 

c 

c 


MAIN PROGRAM 

CHARACTER* 16 FNAME1 , FNAME2 , FNAME3 , FNAME4 
CHARACTER*2 INC1 , INC2 , INC3 , INC4 
OOMMON/BNSN/ EQUIVD.RPM.PEXH.PAIR.TAIR, 

1 CR , BORE , STROKE , CONROD . FPISTE . FPISTA . 

2 EVO , EVC , AVO , AVC , 

3 TWALL , WFUEL , 

4 ACB.ACF, 

5 ANNA , FMEPM , 

6 NPCWF.DELl .DEL2.DEL3, AAN(2) . 

7 VALEXH.CDE.NTEXH, WIDTHE, ALPHEX(50) ,FEXH(50) , 

8 VALAIR . CDA . NTAIR , WIDTHA , ALPAIR(50) . FAIR (50) 
OOMMON/GEN/ALPHAT , ANGEND , ANGRES , AIRPO . EXHPO . IMOLS , 

1 APAIR, APEXH, AREF,AP(2) , 

1 APN(2) , CYCLE .DALPHA, 

1 DREF , DZ , FREF , 

1 IPOWER , IREV , GA , GE , GREF , IUNITL , 

1 IUNITP, IUNITT, IUNITK, IUNITW, IUNITQ, 

1 WAR , PBAR , PBARAB , 

1 PI , PREF . REVENG , REVREF , RPAIR , RPEXH , RP(2) , 

1 RPN(2) ,STEP2,TEXH,TREF, VREF.WREF.XREF, 

1 Z.ZREF, HAIR, HREF, SCMULT.NC 

COMMON/OOMB/ AIRFL , EQUIV . TEXHS , SHPS .FZA.HRD, SCASW , TREFF . SCAEFF . 
1 AIRVMX . EXHVMX , PORTED , STOIC , WOO , BLBACK . FUELF , BHP 

COMMON/OPEN/ FNAME 1 , FNAME2 , FNAME3 , FNAME4 

OOMMON/TEMP/TWEV , TWGH , TWGP , TWGC , TWCR . TWPR , TWCP , TWOC . TWCH , TEXHP , 
1 TWEVG TCI TC2.TC3 

00MM0N/HEAT 1/ QF.QFPC.QC.QCPC, A1 ,A1PC, A2, A2PC, A3.A3PC.A4, A4PC, 

1 QCYL , QCYLPC 

C0MM0N/HEAT2/ QEVG , QEVGPC , QEXHP . QEXHPC , 

1 QTC1 ,QTC1PC,QTC2,QTC2PC,QTC3,QTC3PC,Q0UT,Q0UTPC 
FNAME1=’FILE01.AAA 
FNAME2= ’ FILE01 . BBB 
FNAME3= * FILE01 . OCC 
FNAME4= ’ FILE01 . DDD 
INC1=’02’ 

INC2= ' 03 ’ 

INC3= ’ 04 ’ 

INC4= ’ 05 ’ 

ASSUME COMPRESSOR PRESSURE RATIO AND PRESSURE RATIO ACROSS DIESEL 


P2P 1=6. 0544 
P6P5=1 ./I . 171 


INLET CONDITIONS 
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oooo ooo ooo ooo ooo oo ooo oooo ooo 


Tl=518.7 
Pl=14. 7 

COMPRESSOR OUTLET CONDITIONS 

ETACTT=0 . 85 
P2=P2P1*P1 

DUMMY=P2P1 **0.2857 - 1. 

T2=T1*(DUMMY/ETACTT + 1 . ) 

WC=0 . 24*T1 /ETACTT*DUMMY 

LOSSES TO AFTEROOOLER/AFTERCOOLER INLET CONDITIONS 

P3=0. 98*P2 

P3=P2 

T3=T2 


AFTERCOOLER EXIT CONDITIONS 


E=0.4 

T4=T3-E*(T3-T1) 

P4=0.98*P3 

T4=T3 

P4=P3 


LOSSES TO DIESEL/DIESEL INLET CONDITIONS 


P5=P4 

T5=T4 

PAIR=P5 

TAIR=T5 

PEXH=P6P5*PAIR 

RUN DIESEL PROGRAM 

CALL BENSON 

EXHAUST CONDITIONS 


P6=PEXH 

T6=TEXHS 


LOSSES TO HIGH PRESSURE TURBINE/TURBINE INLET CONDITIONS 

P7=0.98*P6 

T7=T6 

TURBINE OUTLET CONDITIONS - WORK MATCHED TO COMPRESSOR REQUIREMENTS 
ASSUMES 2% LOSS TO BEARING FRICTION 

ICHK=0 


65 



T8=T7 
P8=14. 7 

30 TM=(T8+T7)/2. 

CALL PR0P(EQUIVD , TM , CPX , GX) 

ETAT 1=0.84 

DUMMY=1 . -(P8/P7)**((GX-1 . )/GX) 

T8=T7* ( 1 . -DUMMY*ETAT 1 ) 

WT=ET AT 1 *CPX*T7*DUMMY* ( FUELF+A I RFL ) *3600 . /2544. 
IF ( ICHK . EQ . 1 ) GO TO 40 
ICHK=1 
GO TO 30 
40 CONTINUE 
C 

c 

TWCCF =TWOO 1 . 8-460 . 

TWCPF=TWCP*1 . 8-460 . 

TWCHF=TWCH*1 . 8-460 . 

TWEVF=TWEV*1 . 8-460 . 

TWGHF=TWGH* 1 . 8-460 . 

TWGPF=TWGP* 1 . 8-460 . 

TWPRF=TWPR*1 . 8-460 . 

TWGCF =TWGC* 1 . 8-460 . 

TWCRF=TWCR*1 . 8-460 . 


C 

WRITE(2, 1010) 

1010 FORMAT (/5X, ’ENGINE PART TEMPERATURES’ ,40X, ’HEAT TRANSFER SUMMARY ) 
WRITE(2, 1011) TWGH , TWGHF , QF , QFPC 

1011 FORMAT (/2X, ’GAS SIDE HEAD TEMPERATURE’ .7X.F8. 1 , ’ K’.Fll.l,’ F’, 

1 12X, ’FRICTION GENERATED HEAT TRANSFER’ .3X.F10. 3, FI 1.2, ’ %’) 

WRITE(2. 1012) TWGP , TWGPF , QC , QCPC 

1012 FORMAT (2X, ’GAS SIDE PISTON CROWN TEMPERATURE F7 . 1 , ’ K’.Fll.l, 

1 ’ F’ ,12X, ’CREVICE VOLUME HEAT TRANSFER’ ,F17. 3, FI 1.2, ’ %’) 

WRITE(2, 1013) TWGC, TWGCF. Al.Al PC 

1013 FORMAT (2X, ‘GAS SIDE SLEEVE TEMPERATURE' .5X.F8. 1 , ’ K’.Fll.l,’ F’, 

1 12X, ’HEAT TRANSFER TO CYLINDER WALL’ .F15.3.F11 .2. ’ %’) 

WRITE(2. 1014) TWEV , TWEVF , A2 , A2PC 

1014 F0RMAT(2X, ’EXHAUST VALVE TEMPERATURE 7X , F8 . 1 , ’ K’.Fll.l,’ F’, 

1 12X, 'HEAT TRANSFER TO PISTON’ ,4X,F18.3.F11 .2, ’ %’) 

WRITE(2, 1015) TWPR . TWPRF , A3 , A3PC 

1015 FORMAT (2X. ’PISTON RIM TEMPERATURE’ . 10X.F8. 1 , ’ K’.Fll.l,’ F’. 

1 12X, ’HEAT TRANSFER TO HEAD’ .6X.F18. 3. FI 1.2, ’ %’) 

WRITE(2, 1016) TWCR , TWCRF , A4 , A4PC 

1016 F0RMAT(2X, ’COMPRESSION RING TEMPERATURE ’, 4X . F8 . 1 . ’ K’.Fll.l.' F’, 

1 12X, ’HEAT TRANSFER TO EXHAUST VALVE’ .F15.3.F11 .2. ’ %’) 

WRITE(2, 1017) TWCH , TWCHF , QCYL , QCYLPC 

1017 FORMAT (2X, ’COOLANT SIDE HEAD TEMPERATURE’ .3X.F8. 1 , ’ K’.Fll.l.’ F’, 

1 12X. ’TOTAL HEAT INPUT TO NETWORK* .3X.F15.3.F11 .2, ’ %’) 


WRITE(2, 1018) TWCP.TWCPF 

1018 F0RMAT(2X, 'OIL SIDE PISTON TEMPERATURE ’, 5X , F8 . 1 , ' K’.Fll.l.’ 
WRITE(2. 1019) TWOC, TWOCF 

1019 FORMAT (2X, ’OOOLANT SIDE SLEEVE TEMPERATURE ’. F9 . 1 . ’ K’.Fll.l, 


F*) 

' F’) 
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WRITE(2, 1002) QEVG , QEVGPC , QEXHP , QEXHPC , 

1 QTC1 ,QTC1PC,QTC2,QTC2PC,QTC3,QTC3PC,Q0UT,Q0UTPC 
1002 F0RMAT(69X, ’HEAT TRANSFER TO EXH VALVE GUIDE’ .3X.F10.3.F11 .2, ’ %’/ 

1 69X, ’HEAT TRANSFER TO EXHAUST PORT GASES’ .F10.3.F11 .2, ’ %’/ 

2 69X, 'HEAT TRANSFER TO PISTON COOLING OIL’ .F10.3.F11 .2, ’ %’/ 

3 69X, ‘HEAT TRANSFER TO SLEEVE COOLING OIL’ ,F10.3,F11 .2, ' %’/ 

4 69X. ’HEAT TRANSFER TO HEAD COOLING OIL ’ .F10.3.F11 .2, ’ %’/ 

5 69X, 'TOTAL HEAT TRANSFER FROM NETWORK \F10.3,F11 .2, ’ %’) 

C 

WC=-WO(FUELF+AIRFL)*3600 . /2544 . 

THP=BHP+WC+WT 
BSFC=FUELF/THP*3600 . 

WRITE(2, 1030) BHP.WC.WT.THP.BSFC 
1030 FORMAT (/3X, ’DIESEL BRAKE POWER’ ,F10.2. ’ HP'/ 

1 3X, ’COMPRESSOR POWER’ ,F12.2, ’ HP ’/3X. ’TURBINE POWER’ .F15.2, 

2 ’ HP’/3X. ’TOTAL SHAFT POWER’ ,F11 .2. ’ HP’/3X, ’BSFC’ .20X.F1 1 .4) 

C 

CALL FCLOSE 
STOP 
END 
C 

BLOCK DATA 

COMMON/BNSN/ EQUIVD , RPM . PEXH , PAIR , TAIR , 

1 CR , BORE , STROKE . CONROD , FPISTE . FPISTA . 

2 EVO , EVC , AVO , AVC , 

3 TWALL , WFUEL , 

4 ACB , ACF , 

5 ANNA , FMEPM , 

6 NPCWF , DELI , DEL2 , DEL3 , AAN(2) , 

7 VALEXH.CDE.NTEXH, WIDTHE, ALPHEX(50) ,FEXH(50) , 

8 VALAIR . CDA , NTAIR , WIDTHA , ALPAIR(50> , FAIR(50) 

C 

DATA EQUIVD, RPM, PEXH. PAIR, TAIR /O. 60, 2050. . 134.79, 149.76,896./ 
DATA CR, BORE, STROKE. CONROD, FPISTE, FPISTA /II . 17,6.00,7. 375, 

1 14.75,1. ,1.13/ 

DATA EVO, EVC, AVO. AVC /105. .255. ,130. ,230./ 

DATA TWALL, WFUEL /1460. ,0.0001/ 

DATA ACB, ACF /335.0.10./ 

DATA ANNA. FMEPM /O. 30. 1.25/ 

DATA NPCWF. DELI, DEL2,DEL3,AAN(1).AAN(2) /-I ,0.01457.0.667,0.4. 

1 0.4,0. 5/ 

DATA VALEXH.CDE, WIDTHE /O. ,0.85,0.20/ 

DATA VALAIR, CDA, WIDTHA /O. ,0.85,0.36/ 

END 

C 

SUBROUTINE PROP(EQUIV,T,CPX,GX) 

IF(EQUIV.GT.O.) GO TO 10 
X1=0. 

X2=0. 

X3=0. 21 
X4=0.79 
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GO TO 20 

10 CARBON =0.85561 

HC= 1 2 . 0 1 1 / 1 . 008* ( 1 . -CARBON ) /CARBON 
AST0=(2. +HC/2 . )/2 . /O . 21 
Yl=l . 

Y2=HC/2 

Y3=AST0/EQUI VO . 2 1 - 1. - HC/4. 

Y4=ASTT0/EQUIV * 0.79 

YT0T=Y1+Y2+Y3+Y4 

X1=Y1/YT0T 

X2=Y2/YT0T 

X3=Y3/YT0T 

X4=Y4/YTOT 

20 CP1=16. 2-6.53E3/T+1 .41E6/T**2 
CP2= 1 9 . 86-597 . /SQRT (T)+7500./T 
CP3=11. 515-172. /SQRT(T)+1530./T 
CP4=9 . 47-3 . 47E3/T+ 1 . 1 6E6/T**2 
CPX=X 1 OP 1 +X20P2+X30P3+X40P4 
AVM=X 1 *44 . 0 1 +X2* 1 8 . 0 1 6+X3*32 . 00+X4*2S . 0 1 
CPX=CPX/AVM 
RX= 1 545 . /778 . /AVM 
CVX=CPX-RX 
GX=CPX/CVX 
RETURN 
END 
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10.3 Garrett CCE program listing 
MAIN PROGRAM 

CHARACTER* 16 FNAME 1 , FNAME2 , FNAME3 , FNAME4 
CHARACTER*2 INC1 . INC2 , INC3 . INC4 
OOMMON/BNSN/ EQUIVD.RPM.PEXH.PAIR.TAIR, 

1 CR . BORE . STROKE , CONROD . FPISTE . FPISTA . 

2 EVO.EVC.AVO.AVC, 

3 TWALL.WFUEL. 

4 ACB.ACF. 

5 ANNA , FMEPM . 

6 NPCWF.DELl .DEL2.DEL3, AAN(2) , 

7 VALEXH.CDE.NTEXH.WIDTHE, ALPHEX(50) ,FEXH(50) , 

8 VALAIR , CDA , NTAIR , WIDTHA , ALPAIR(50) , FAIR(50) 
COMMON/GEN/ALPHAT , ANGEND , ANGRES , AIRPO , EXHPO . IMOLS . 

1 APAIR, APEXH, AREF, AP(2) , 

1 APN(2), CYCLE, DALPHA, 

1 DREF.DZ.FREF, 

1 I POWER, IREV ,GA,GE, GREF, IUNITL. 

1 IUNITP, IUNITT, IUNITK, IUNITW, IUNITQ, 

1 WAR , PBAR , PBARAB , 

1 PI , PREF , REVENG , REVREF , RPAIR , RPEXH ,RP(2) , 

1 RPN (2) , STEP2 , TEXH , TREF , VREF , WREF , XREF , 

1 Z.ZREF, HAIR, HREF, SCMULT.NC 

COMMON/COMB/ AIRFL , EQUIV , TEXHS , SHPS .FZA.HRD, SCASW , TREFF , SCAEFF , 
1 AIRVMX , EXHVMX , PORTED , STOIC , WCO , BLBACK , FUELF , BHP 
COMMON/OPEN/ FNAME 1 , FNAME2 , FNAME3 , FNAME4 

OOMMON/TEMP/TWEV , TWGH , TWGP , TWGC , TWCR , TWPR , TWCP , TWOC , TWCH , TEXHP , 
1 TWEVG TC 1 TC2 TC3 

COMMON/HEAT 1 / QF.QFPC.QC.QCPC, A1 , A1PC, A2, A2PC, A3, A3PC, A4, A4PC, 

1 QCYL , QCYLPC 

C0MM0N/HEAT2/ QEVG , QEVGPC , QEXHP , QEXHPC , 

1 QTC1 , QTC1PC , QTC2 . QTC2PC , QTC3 . QTC3PC , QOUT . QOUTPC 

FNAME1=’FILE01 .AAA 
FNAME2=’FILE01 .BBB 
FNAME3= ’ FILE01 . OOC 
FNAME4- ’ FILEO 1 . DDD 
INC1=’02’ 

INC2= ' 03 ’ 

INC3= ’ 04 ’ 

INC4= ’ 05 ' 

C 

C ASSUME COMPRESSOR PRESSURE RATIO AND PRESSURE RATIO ACROSS DIESEL 
C 

P2P1=10.61 
P6P5=1 ./I . Ill 
C 

C INLET CONDITIONS 
C 


Tl=518.7 



non ooo ooo ooo ooo ooo oooo ooo 


Pl=14.7 


COMPRESSOR OUTLET CONDITIONS 

ETACTT=0.79 

P2=P2P1*P1 

DUMMY=P2P1**0. 2857 - 1. 

T2=T 1 * ( DUMMY/ET ACTT + 1 . ) 

WC=0 . 24*T 1 /ETACTT*DUMMY 

LOSSES TO AFTEROOOLER/AFTEROOOLER INLET CONDITIONS 

P3=0.98*P2 

P3=P2 

T3=T2 


AFTEROOOLER EXIT CONDITIONS 


E=0. 4 

T 4=T3-E* ( T3-T 1 ) 

P4=0.98*P3 

LOSSES TO DIESEL/DIESEL INLET CONDITIONS 

P5=P4*0 . 98 

T5=T4 

PAIR=P5 

TAIR=T5 

PEXH=P6P5«PAIR 

RUN DIESEL PROGRAM 

CALL BENSON 

EXHAUST CONDITIONS 


P6=PEXH 

T6=TEXHS 


LOSSES TO HIGH PRESSURE TURBINE/TURBINE INLET CONDITIONS 

P7=0.98*P6 

T7=T6 


AFTERBURNER OUTLET CONDITIONS 


F=0.0 

TIN=T6 

CALL TBURN(EQUIVD,F, TIN, TOUT) 
T7=T0UT 

XFUEL=F*(AIRFL+FUELF) 
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ooo ooo o o o o 


TURBINE OUTLET CONDITIONS - WORK MATCHED TO COMPRESSOR REQUIREMENTS 
ASSUMES 2% LOSS TO BEARING FRICTION 


ICHK=0 

T8=T7 

30 TM=(T8+T7)/2. 

CALL PROP(EQUIVD , TM , CPX , GX) 

ET AT 1=0. 845 
WT1=1.02*WC 

P8=P7*K1 ■-WT1/(ETAT1*TX*T7))**(GX/(GX-1 . )) 

DUMMY=1 .-(P8/P7)**( (GX-1 . )/GX) 

T8=T7*( 1 . -DUMMYwETATl ) 

IF(ICHK.EQ.l) GO TO 40 
ICHK=1 
GO TO 30 
40 CONTINUE 

LOSSES BETWEEN TURBINES/POWER TURBINE INLET CONDITIONS 

P9=0 . 98*P8 
T9=T8 


POWER TURBINE OUTLET CONDITIONS 

ICHK=0 
T10=T9 

50 TM=(T9+T10)/2. 

CALL PROP ( EQU I VD , TM , CPX , GX ) 

WRITE(6,*) ’ CPX , GX ’ , CPX , GX 
P10=14.7 
ETAT2=0 . 848 

DUMMY=1 .-(P10/P9)**((GX-1 . )/GX) 

T10=T9*( 1 . - DUMMY *ETAT2) 

IF(ICHK.EQ.l) GO TO 60 
ICHK=1 
GO TO 50 
60 CONTINUE 

WT2=ETAT2*CPX*T9*DUMMY * ( FUELF+ A IRFL ) *3600 . /2544 . 
C 

WRITE(6,*) ’T9,P9,T10,P10’ .T9.P9.T10. P10 
C 

TW0CF=TW0Ol . 8-460 . 

TWCPF=TWCP*1 . 8-460 . 

TWCHF=TWCH*1 . 8-460 . 

TWEVF=TWEV*1 . 8-460 . 

TWGHF=TWGHxl . 8-460 . 

TWGPF=TWGP*1 . 8-460 . 

TWPRF=TWPR« 1 . 8-460 . 

TWGCF=TWG01 . 8-460 . 

TWCRF =TWCR* 1 . 8-460 . 
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c 

WRITE(2, 1010) 

1010 FORMAT (/5X, ’ENGINE PART TEMPERATURES ’ , 40X . ' HEAT TRANSFER SUMMARY’) 
WRITE(2, 1011) TWGH . TWGHF , QF . QFPC 

1011 FORMAT(/2X, ’GAS SIDE HEAD TEMPERATURE’ .7X.F8. 1 . ’ K’.Fll.l.’ F’. 

1 12X, ’FRICTION GENERATED HEAT TRANSFER’ .3X.F10. 3. Fll .2, ’ %’) 

WRITE(2, 1012) TWGP , TWGPF , QC , QCPC 

1012 FORMAT (2X, ’GAS SIDE PISTON CROWN TEMPERATURE ’. F7 . 1 , ’ K’.Fll.l, 

1 ’ F* .12X. ’CREVICE VOLUME HEAT TRANSFER’ ,F17. 3, Fll .2. ’ %’} 

WRITE(2, 1013) TWGC.TWGCF, A1 , A1PC 

1013 FORMAT (2X, ’GAS SIDE SLEEVE TEMPERATURE’ .5X.F8. 1 . ’ K’.Fll.l.’ F’. 

1 12X, ’HEAT TRANSFER TO CYLINDER WALL’ .F15.3.F11 .2. ’ %’) 

WRITE(2, 1014) TWEV . TWEVF , A2 , A2PC 

1014 FORMAT (2X. ’EXHAUST VALVE TEMPERATURE ’. 7X . F8 . 1 . ’ K’.Fll.l,’ F’, 

1 12X, 'HEAT TRANSFER TO PISTON’ .4X.F18.3.F11.2, ' %’) 

WRITE(2, 1015) TWPR , TWPRF , A3 , A3PC 

1015 F0RMAT(2X, ’PISTON RIM TEMPERATURE’ , 10X.F8. 1 , ’ K’.Fll.l,’ F’, 

1 12X. ’HEAT TRANSFER TO HEAD’ .6X.F18. 3. Fll. 2. ’ %’) 

WRITE(2, 1016) TWCR , TWCRF , A4 , A4PC 

1016 FORMAT (2X, 'COMPRESSION RING TEMPERATURE ’, 4X . F8 . 1 . ’ K’.Fll.l,’ F', 

1 12X, ’HEAT TRANSFER TO EXHAUST VALVE’ .F15.3.F11 .2, ’ %’) 

WRITE(2, 1017) TWCH , TWCHF , QCYL , QCYLPC 

1017 FORMAT (2X, ’COOLANT SIDE HEAD TEMPERATURE ’. 3X . F8 . 1 , ’ K’.Fll.l,’ F’, 

1 12X, ’TOTAL HEAT INPUT TO NETWORK' .3X.F15.3.F11 .2, ’ %’) 

WRITE(2, 1018) TWCP.TWCPF 

1018 F0RMAT(2X, ’OIL SIDE PISTON TEMPERATURE ’, 5X . F8 . 1 . ’ K’.Fll.l,’ F’) 
WRITE(2, 1019) TWCC.TWOCF 

1019 FORMAT ( 2X, ’OOOLANT SIDE SLEEVE TEMPERATURE ’, F9 . 1 , ’ K’.Fll.l,’ F’) 
WRITE(2, 1002) QEVG , QEVGPC , QEXHP , QEXHPC , 

1 QTC1,QTC1PC.QTC2.QTC2PC,QTC3,QTC3PC,Q0UT,Q0UTPC 
1002 FORMAT (69X, ’HEAT TRANSFER TO EXH VALVE GUIDE’ ,3X,F10.3,F11 .2, ’ %’/ 

1 69X, 'HEAT TRANSFER TO EXHAUST PORT GASES’ .F10.3.F11 .2, ’ %’/ 

2 69X, ’HEAT TRANSFER TO PISTON OOOLING OIL’ .F10.3.F11 .2, ’ %'/ 

3 69X, ’HEAT TRANSFER TO SLEEVE OOOLING OIL* .F10.3.F11 .2, ’ %’/ 

4 69X, ’HEAT TRANSFER TO HEAD OOOLING OIL ’ .F10.3.F11 .2, ’ %’/ 

5 69X, ’TOTAL HEAT TRANSFER FROM NETWORK ’ .F10.3.F11 .2, ’ %’) 

C 

FANP=51./6. 

THP=BHP+WT2-FANP 

BSFC= ( FUELF+XFUEL ) /THP*3600 . 

WRITE(2, 1030) BHP , WT2 , FANP , THP . BSFC 
1030 FORMAT (/3X, ’DIESEL BRAKE POWER’ ,F10.2, ’ HP’/ 

1 3X. ’TURBINE POWER’ ,F15. 2. ’ HP’/ 

2 3X, ’FAN POWER ’.F15.2. 

3 ’ HP’/3X. ’TOTAL SHAFT POWER’ , Fll .2, ’ HP’/3X, ’BSFC’ , 15X.F11 .4) 
WRITE(6,*) 'BSFC = ' .BSFC 

C 

CALL FCLOSE 
STOP 
END 
C 


72 



BLOCK DATA 

OOMMON/BNSN/ EQUIVD . RPM , PEXH , PAIR , TAIR , 

1 CR , BORE , STROKE , CONROD , FPISTE . FPISTA , 

2 EVO . EVC , AVO , AVC , 

3 TWALL.WFUEL, 

4 ACB.ACF, 

5 ANNA , FMEPM , 

6 NPCWF , DELI , DEL2 , DEL3 , AAN(2) . 

7 VALEXH,CDE,NTEXH, WIDTHE, ALPHEX(50) ,FEXH(50) , 

8 VALAIR , CDA , NTAIR , WIDTHA , ALPAIR(50) . FAIR(50) 
Cxxxxxxxxxxxxxxxxx 

DATA EQUIVD. RPM, PEXH, PAIR, TAIR /0. 68, 6122. . 134.79, 149.76,896./ 
DATA CR, BORE, STROKE, OONROD, FPISTE, FPISTA /9. 1712,3. 101 ,2.940, 

1 7.168,1. ,1.13/ 

DATA EVO, EVC. AVO, AVC /90. ,239. , 126. .234./ 

DATA TWALL.WFUEL / 1460. ,0.000145332/ 

DATA ACB.ACF /345..10.5/ 

DATA ANNA , FMEPM /O . 36 , 1 . 0/ 

DATA NPCWF, DELI , DEL2 , DEL3 , AAN ( 1 ) ,AAN(2) /-I .0.01457,0.667,0.4, 

1 0.4, 0.5/ 

DATA VALEXH , CDF , NTEXH /10. ,0.8864, 13/ 

DATA VALAIR, CDA. WIDTHA /O. ,0.8,0.5542/ 

DATA ALPHEX/O . 00 , 10.97,25.00,31.61,41.00.52.25,72.89,90.95, 104.00, 
1 111.59,119.00,132.23, 149.00,37*0./ 

DATA FEXH/O. 0,0. 00093264, 0.00364, 0.0059 1222, 0.00899. 0.0 106201 , 

1 0.01175263,0.01111131 .0.00923,0.00717895.0.00461 ,0.00181743,0. , 

2 37*0./ 

END 

C 

SUBROUTINE PROP(EQUIV,T,CPX,GX) 

IF(EQUIV.GT.O. ) GO TO 10 
XI =0. 

X2=0. 

X3=0. 21 
X4=0.79 
GO TO 20 

10 CARBON=0. 85561 

HC= 1 2 . 0 1 1 / 1 . 008* ( 1 . -CARBON ) /CARBON 

AST0=(2. +HC/2 . )/2 . /O . 21 

Yl=l. 

Y2=HC/2 . 

Y3=ASTT0/EQUIV*0.21 - 1. - HC/4. 

Y 4=AST0/EQU IV * 0.79 

YT0T=Y1+Y2+Y3+Y4 

XUY1/YT0T 

X2=Y2/YT0T 

X3=Y3/YT0T 

X4=Y4/YT0T 

20 CPU16.2-6.53E3/T+1 .41E6/T**2 

CP2= 19.86-597. /SQRT ( T ) +7500 . /T 
CP3=11. 515-172. /SQRT(T)+1530./T 



CP4=9 . 47-3 . 47E3/T+ 1 . 16E6/T**2 

CPX-X 1 *CP 1 +X2*CP2+X3*CP3+X4*CP4 

AVM=X 1 *44 . 0 1 +X2* 18.01 6+X3*32 . 00+X4*28 . 0 1 

CPX=CPX/AVM 

RX=1545 . /778 . /AVM 

CVX=CPX-RX 

GX=CPX/CVX 

RETURN 

END 

SUBROUTINE TBURN(EQUIV.F.T1 ,T2) 

HV__ 18600 

CALL PR0P(EQUIV,T1 ,CP1 ,G1 ) 

T537=537 

CALL PROP(EQUIV , T537 , CPO , GO) 

CPAVE= (CP 1 +CPO ) /2 . 

DELH 1 =CPAVE* ( T 1 -T537 ) 

DELH=DELH1+HV*F 

CARBON =0.85561 

STOI C= 1 ./(( CARBON/ 1 2 . 0 1 + ( 1 . -CARBON )/( 1 . 008*4 . )) * ( 32 . 0+3 . 76*28 . 0 ) ) 
EQUIV2=(F+1.)*EQUIV + F/STOIC 
WRITE(6,*) ’ STOIC. EQUIV2’ . STOIC , EQUIV2 
T0LD=2300 . 

CALL PR0P(EQUIV2 , TOLD , CP2 . G2) 

T537=537. 

CALL PR0P(EQUIV2 . T537 , CPO .GO) 

CPAVE= (CP2+CP0 ) /2 . 

DELH2=CPAVE*(T0LD-T537)*( l . +F) 

Z0LD=DELH2-DELH 
TNEW=TOLD* 1 . 05 

CALL PROP ( EQU IV2 . TNEW . CP2 . G2 ) 

CPAVE= ( CP2+CP0 ) /2 . 

DELH2=CPAVE« ( TNEW-T537 ) ^ ( 1 .+F) 

ZNEW=DELH2-DELH 

DELT=ZNEW« ( TNEW-TOLD ) / ( ZOLD-ZNEW ) 

T2=TNEW+DELT 

WRITE(6.«) ’ ZNEW . T2 ’ , ZNEW , T2 

IF(ABS(DELT).LT. 0.001) GO TO 50 

T0LD=TNEW 

TNEW=T2 

ZOLD=ZNEW 

GO TO 25 

WRITE(6.«) ’ DELH. DELH2’ .DELH. DELH2 

RETURN 

END 


10 . 4 Benson Program Listing 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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VARIABLE NAMES : 

ACBB = DUMMY VARIABLE USED TO STORE ACB, START OF COMBUSTION 
ACFF = DUMMY VARIABLE USED TO STORE ACF , END OF COMBUSTION 
AIR = AN INDICATOR OF WHETHER INTAKE VALVE/PORT IS OPEN OR 
CLOSED DURING THE CYCLE (=0 WHEN CLOSED, =10 WHEN 
OPEN. 

ALPAIR = ARRAY OF CRANK POSITIONS CORRESPONDING TO THE INTAKE 
VALVE FLOW AREAS GIVEN IN FAIR, DEGREES (INPUT) 

ALPHAT = CUMMULATIVE TOTAL CRANKANGLE DURING INTEGRATION 
ALPHEX = ARRAY OF CRANK POSITIONS CORRESPONDING TO THE EXHAUST 
VALVE FLOW AREAS GIVEN IN FEXH, DEGREES (INPUT) 

ANGEND = MAXIMUM NUMBER OF CRANK DEGREES ALLOWED BEFORE 
INTEGRATION IS TERMINATED 

AAN ( 1) = FUEL INJECTION PARAMETER FOR WHITEHOUSE-WAY MODEL 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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AAN ( 2 ) = FUEL INJECTION PARAMETER FOR WHITEHOUSE-WAY MODEL 
ACB = BEGINNING OF INJECTION WHEN NPCWF . LT. 0 

BEGINNING OF HEAT RELEASE WHEN NPCWF. GT. 0 
ACF = END OF INJECTION WHEN NPCWF. LT. 0 

END OF HEAT RELEASE WHEN NPCWF . GT . 0 
AVC = INTAKE VALVE/PORT OPENING ANGLE, DEGREES 
AVO = INTAKE VALVE/PORT CLOSING ANGLE, DEGREES 
AREF = REFERENCE VALUE FOR SPEED OF SOUND (BASED ON TREF) 
ANNA = PARAMETER FOR ANNAND OR WOSCHNI HEAT TRANSFER MODEL 
ANNB = PARAMETER FOR ANNAND OR WOSCHNI HEAT TRANSFER MODEL 
ANNC = PARAMETER FOR ANNAND OR WOSCHNI HEAT TRANSFER MODEL 
BETA = PARAMETER IN WATSON COMBUSTION MODEL 
CARBON = MASS FRACTION OF CARBON IN FUEL 
CD1 = PARAMETER IN WATSON COMBUSTION MODEL 
CD2 = PARAMETER IN WATSON COMBUSTION MODEL 
CYCLE = 2 FOR TWO CYCLE ENGINE, = 4 FOR FOUR CYCLE ENGINE 
COEFFA = COEFFICIENT FOR POLYNOMIAL FIT TO JANAF TABLES 
COEFFB = COEFFICIENT FOR POLYNOMIAL FIT TO JANAF TABLES 
COEFFC = COEFFICIENT FOR POLYNOMIAL FIT TO JANAF TABLES 
GOEFFD = COEFFICIENT FOR POLYNOMIAL FIT TO JANAF TABLES 
CPI = PARAMETER IN WATSON OCMBUSTION MODEL 
CP2 = PARAMETER IN WATSON COMBUSTION MODEL 
CYCLE = 2 FOR TOO STROKE ENGINE, = 4 FOR A POUR STROKE ENGINE 
SURFC = SURFACE AREA OF CYLINDER WALL EXPOSED TO COMBUSTION 
DELI = PARAMETER FOR EITHER WIEBE OR WHITEHOUSE-WAY MODEL 
DEL2 = PARAMETER FOR EITHER WIEBE OR WHITEHOUSE-WAY MODEL 
DEL3 = PARAMETER FOR EITHER WIEBE OR WHITEHOUSE-WAY MODEL 
CHEAT = HEAT TRANSFER RATE FROM CYLINDER CURING GAS EXCHANGE 
DQF = HEAT ADDITION TO CYLINDER DUE TO COMBUSTION 
DQCR = HEAT TRANSFER FROM CYLINDER DUE TO RADIATION 
DQCY = HEAT TRANSFER FROM CYLINDER DUE TO CONVECTION 
ETI = INTERNAL ENERGY OF GASES IN CYLINDER AT START OF 
TIME STEP 

ETF = INTERNAL ENERGY OF GASES IN CYLINDER AT END OF 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


TIME STEP 

EVC = EXHAUST VALVE/PORT OPENING ANGIE , DEGREES 
EVO = EXHAUST VALVE/PORT CLOSING ANGIE, DEGREES 
EXH = AN INDICATOR OF WHETHER THE EXHAUST VALVE/PORT IS OPEN 
OR CLOSED DURING THE CYCLE. (=0 WHEN CIOSED, =10 WHEN 
OPEN) 

FCYL = CROSSECnONAL AREA OF CYLINDER 
PMEP1 = PARAMETER IN FMEP CORRELATION 
FMEP2 = PARAMETER IN FMEP CORRELATION 
FMEP3 = PARAMETER IN FMEP CORRELATION 

FPISTA = RATIO OF SURFACE AREA OF PISTON TOP TO CROSSECTI ONAL 
AREA OF CYLINDER 

FPISTE = RATIO OF SURFACE AREA OF HEAD EXPOSED TO COMBUSTION 
CHAMBER TO CROSSECTIONAL AREA OF CYLINDER 
FAIR = ARRAY SPECIFYING THE GEOMETRIC INTAKE VALVE FLOW AREA 
FEXH = ARRAY SPECIFYING THE GEOMETRIC EXHAUST VALVE FLOW AREA 
FREF = A REFERENCE AREA USED TO NORMALIZE AREAS IN THE PROGRAM 
FZA = FUEL/AIR MASS RATIO 
GA = RATIO OF SPECIFIC HEATS FOR AIR 
GE = RATIO OF SPECIFIC HEATS FOR EXHAUST PRODUCTS 
GREF = REFERENCE RATIO OF SPECIFIC HEATS (INFUT) 

HEATG = TOTAL HEAT TRANSFER DURING GAS EXCHANGE 
NPCWF = PARAMETER WHICH DETERMINES WHICH OOMBUSTION MODEL IS 
USED IN THE PROGRAM 

WHEN NPCWF . LT . 0 WHITEHOUSE-WAY MODEL 
WHEN NPCWF. EQ.O AND DEL3.EQ.0 AVL MODEL 
WHEN NPCWF. EQ.O AND DEL3.GT.0 WIEBE MODEL 
WHEN NPCWF . GT . 0 WATSON MODEL 

NTAIR = NUMBER OF POINTS PROVIDED IN ALPAIR AND FAIR TO SPECIFY 
THE INTAKE VALVE FIOW AREA. FAIR (NTAIR) , ALPAIR (NTAIR) 
NTEXH = NUMBER OF POINTS PROVIDED IN ALFHEX AND FEXH TO SPECIFY 
THE EXHAUST VALVE FLOW AREA. FEXH (NTEXH) , ALFHEX (NTEXH) 

P02 = PARTIAL PRESSURE OF OXYGEN 

PORTED = 1 FOR LOOP SCAVENGED, = 2 FOR UNIFICW SCAVENGED 
STOIC = STOICHIOMETRIC FUEL/ AIR RATIO 
SURFP = SURFACE AREA OF PISTON 

SURFH = SURFACE AREA OF HEAD EXPOSED TO OOMBUSTION CHAMBER 
TAIR = TEMPERATURE OF AIR ENTERING CYLINDER THROUGH INTAKE 
TEXH = EXHAUST TEMPERATURE (INITIAL ESTIMATE SPECIFIED AS INPUT) 
TREF = REFERENCE TEMPERATURE (INPUT) 

TWALL = CYLINDER WALL TEMPERATURE 

VALAIR = 0 IF ENGINE HAS INTAKE PORTS, = 10 IF ENGINE HAS 
INTAKE VALVES 

VAIEXH = 0 IF ENGINE HAS EXHAUST PORTS, = 10 IF ENGINE HAS 
EXHAUST VALVES 

WABU = MASS OF AIR REQUIRED TO BURN FUEL UNDER STOICHIOMETRIC 
CONDITIONS, MASS UNITS 

WCA = AIR RETAINED IN CYLINDER PER CYCLE, MASS UNITS 
WCHG = MASS OF AIR PLUS RESIDUAL IN CYLINDER PER CYCLE, MASS UNITS 
WOO = TOTAL MASS FLOW RATE THROUGH EXHAUST VALVE, MASS/TIME 
WIDIHA = FRACTION OF CYLINDER BORE OCCUPIED BY INTAKE PORTS 
WIDIHE = FRACTION OF CYLINDER BORE OCCUPIED BY EXHAUST PORTS 
WFUEL = MASS OF FUEL INJECTED PER CYCLE, MASS UNITS 
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C WMW = MOLECULAR WEIGHT 
C WN(I) = NUMBER OF MOLES OF SPECIES I 
C WNA = NUMBER OF MOLES OF AIR TRAPPED IN CYLINDER 
C WNI (I) = NUMBER OF MOLES OF EACH SPECIES PRESENT AT START OF 
C TIME STEP 

C WNF(I) = NUMBER OF MOLES OF EACH SPECIES PRESENT AT END OF 
C TIME STEP 

C WNP = INCREASE IN NUMBER OF MOLES IN CYLINDER DUE TO 
C COMBUSTION (MOLS AIR —MOLES PRODUCTS) 

C WNR = NUMBER OF MOLES OF RESIDUAL TRAPPED IN CYLINDER 

C WNTF = TOTAL NUMBER OF MOLES OF GASES IN CYLINDER, GMOLES 

C WORKGA = WORK DONE WHEN VALVES ARE CLOSED 
C WORKGE = WORK DONE DURING GAS EXCHANGE 
C WRES = MASS OF RESIDUAL GASES IN CYLINDER, MASS UNITS 
C WTOT = TOTAL MASS IN CYLINDER, AIR+RESIDUAL4-FUEL, MASS UNITS 
C 
C 

c 

SUBROUTINE BENSON 

CYLINDER PRESSURE, TEMPERATURE, VOLUME ON DIESELS . OUT 
IMPLICIT REAL ( A-H, O-Z ) , INTEGER ( I— N) 

COMMON/BNSN/ EQUIVZ , RPZ , PEXZ , PAIZ , TAIZ , 

1 CZ , BORZ , STROKZ , OONROZ , FPISTY , FPISTZ , 

2 EVOZ , EVCZ , AVOZ , AVCZ , 

3 TWALLZ,WFUEZ, 

4 ACBZ , ACFZ , 

5 ANNAZ , FMEPMZ , 

6 NPCWFZ,DEL1Z,DEL2Z,DEL3Z,AANZ(2) , 

7 VALEXZ,CDEZ,NTEXHZ,WIDTZE,ALFHEZ(50) ,FEXHZ(50) , 

8 VALAIZ, CDAZ ,NTAIRZ, WIDTZA, ALPAIZ (50) , FAIRZ (50) 
OQMMON/XXX/ EQUTVD, RIM, PEXH, PAIR, TAIR, 

1 CR , BORE , STROKE , CONROD , FPISTE , FPISTA , 

2 EVO,EVC, AVO,AVC, 

3 TWALL,WFUEL, 

4 ACB,ACF, 

5 ANNA, FMEEW, 

6 NPCWF,DEL1,DEL2,DEL3,AAN(2) , 

7 VALEXH,CDE,NTEXH,WIDrHE,ALFHEX(50) ,FEXH(50) , 

8 VALAIR,CDA,NTAIR,WIDrHA, ALPAIR(50) , FAIR (50) 

COMMON /GEN/ALPHAT , ANGEND , ANGRES , AIRPO , EXHPO , IMOLS , 

1 APATR,APEXH,AREF, AP(2) , 

1 APN(2) , CYCLE, DALFHA, 

1 DREF,DZ,FREF, 

1 IPOWER, XREV,GA,GE,GREF,IUNITL, 

1 IUNITP, IUNTTT, IUNITK, IUNITW, IUNITQ, 

1 WAR,PBARAB, 

1 PI,PREF,REVENG,REVREF,RPAIR,RPEXH,RP(2) , 

1 RPN(2) ,STEP2,TEXH,TREF,VREF,WREF,XREF, 

1 Z,ZREF, HAIR, HREF, SCMULT,NC 

OOMMON/CYI/ACSB , ALFHA , ALFHAE , ANNB , ANNC , AC , 

1 ACN,ACR,ALPHAC, 

1 CALVAL, CARBON, 

1 COEFFA(4) ,00EFFB(4) ,COEFFC(4) ,00EFFD(4) , 

1 COEFFE ( 4 ) , COEFFZ ( 4 ) , CRANK, DCYL, CHEAT, 
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1 DRC,WCTL,DWCIN,DWCP(2) ,EWCYL, 

1 FCYL,HEATG, 

1 RJRITY, PVEL,PCR, PCYLT, 

1 PORTS, RC,RCN,RCR,SURFC, 

1 TCR,TCYLT, VSWEPT, 

1 VCYL,WCAIR / WCIN, 

1 WCOUT,WCYL,VMW(4) ,WPaTT(4) ,WCYIR, 

1 WN (4) , WORKGA , WORKGE , XCRA , XSTA , 

1 NCD , WOTTJT , WCINT , 

1 TPIST,THEAD,TEXV, CID 

aa^MON/CXXffi/AIRFL , EQCJIV , TEXHS , SHPS , FZA , HRD , SCASW , TREFF , SCAEFF , 
+ AIRVMX , EXHVMX , PORTED , STOI C , WOO , BLBACK , FUELF , BHP 
OCMMON/HEAT/ NHEAT , QSLEV1 , QPIST1 , QHEAD1 , QEXHV1 , QTOT1 , 

1 QSLEV2 , QPIST2 , QHEAD2 , QEXHV2 , QTOT2 , 

2 HC1,HC2,HC3,HC4,HC5,HC6,HC7 ,HC8, 

3 HCTl , HCT2 , HCT3 , HCT4 , HCT5 , HCT6 , HCT7 , HCT8 
COMMON/STOP/ ISTOP 

WRITE(6, *) 'INSIDE BENSON' 

EQUIVD=EQUTVZ 

RPW=RPZ 

PEXH=PEXZ 

PAIR=PAIZ 

TAIR=TAIZ 

CR=CZ 

BORE=BORZ 

STROKE=STROKZ 

OONROD=OONROZ 

FPISTE=FPISTY 

FPISTA=FPISTZ 

EVO=EVOZ 

EVC=EVCZ 

AVO=AVOZ 

AVC=AVCZ 

IWALI^IWALLZ 

WFUED=OTUEZ 

ACB=ACBZ 

ACF=ACFZ 

ANNA=ANNAZ 

PMEFM=FMEFMZ 

NPCWF=NPCWFZ 

DEL1=DEL1Z 

DEL2=DEL2Z 

DEL3=DEL3Z 

AAN (1) =AANZ (1) 

AAN ( 2 ) =AANZ ( 2 ) 

VALEXH=VALEXZ 

CDE=CDEZ 

NTEXH=NTEXHZ 

WIDIHE=WIDTZE 

VAIAIR=VAIAIZ 

CDA=CDAZ 

NTAIR=NrAIRZ 

WIDIHA=WIDTZA 

DO 444 1=1,50 
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ALPHEX (I) =ALFHEZ (I) 

FEXH(I)=FEXHZ (I) 

ALPAIR(I)=ALPAIZ(I) 

444 FAIR ( I ) =FAIRZ ( I ) 

C 

C 

ISTOP=0 
ZERO=0 • 0 
0HE=1.0 
PI=3. 1415927 
C 

C SELECT HEAT TRANSFER MODEL 
C 

NHEAT=1 

C 

C OPEN OUTPUT DATA FILES 
C 

CALL FOFEN 

WRITE (6, *) 'RETURNED FROM EOFEN' 
c 

c This section gets and prints the time and date of this run on the outputs 
c 

C CALL GEITIM (IHOUR, 3MINUT, IDUM1 , IDUM2 ) 

C CALLGETDAT (IYEAR, IMONTH,IDAY) 

C DO 2 IZZZ=2 , 4 

C WRITE (IZZZ,l)IMONTH,IDAY, IYEAR, IHOUR, IMINUT 

C 1 FORMAT (' DATE: ' ,12, '/' ,12, '/' ,I4,15X, 'TIME: ',12,':', 12) 

C 2 CONTINUE 
C 

TEXHS = 0. 

SHPS = 0. 

FZA = 0. 

EQUTV = .5 
AIRVMX = 0. 

EXHVMX = 0. 

C 

CALL INPUT 

WRITE (6,*) 'RETURNED FROM INPUT' 

C 

CALL INEUT2 

WRITE (6,*) 'RETURNED FROM INPUT2 ' 

C 

CALL HEADER 

WRITE (6,*) 'REIURNED FROM HEADER' 
c 

WRITE (3, 2048) 

2048 FORMAT (4X, 5HANGLE,5X, 8HPRESSURE, 3X, 11HTEMPERATORE, 4X, 

+ 5HGAMMA, 5X, 6HVOLUME, 5X, 9HMEAN TEMP/) 

JJJJ=1 

Z=0.000 

CALL UNITP ( IUNITP , 1 , PREF , PBARAB) 

C SET UP TEMPERATURES TO S. I. UNITS 
CALL UNITT ( IUNITT , 1 , TREF) 

AREF=SQRT (287 . 1*GREF*TREF) 
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C ENTER CALCULATION FIRST TIME 
XREF=1 . 

EREF=1. 

FREF=1. 

VREF=1. 

WREF=GREF*FREF*VREF* 1 . 0E5/AREF**2 
ZREF=360 . 0*XREF*REVREF/AREF 
REVENG=REVENG/REVREF 
CALL UNITT ( IUNITT , 1 , TEXH) 

CALL UNFIT ( IUNITT , 1 , TAIR) 

CALL UNITP ( IUNITP , 1 , FEXH , FBARAB) 

CALL UNITP ( IUNITP, 1, PAIR, PBARAB) 

CALL BAL(0. ,0. ,0. ,0. , TAIR, GA, HAIR) 

CALLBAL(0. ,0. ,0. ,0. ,TREF, CUM, HREF) 

WRITE(6,*) 'ENTERING LOOP' 

C 

6 RPEXH=PEXH/FREF 

rpair=pair/pref 

APEXH= (SQRT (0 . 287E3*GE*TEXH) ) /AREF 
AAEXH=APEXH/RPEXH* * ( (GEHL . 0) / (2 . 0*GE) ) 

APAIR= (SQRT ( 0 . 287E3*GA*TAIR) ) /AREF 
C 

DZ = STEP2/ (ZREF*REVENG) 

Z = Z+DZ 

ALFHAT = Z * ZREF*REVENG 
C 

IF (JJJJ.EQ. 1) CALL INITLZ 
C 

CALL CYLNDR(JJJJ) 

CHECK TO SEE IF CALCULATION IS COMPUTED 
IF(ISTOP.EQ.l) GO TO 999 

CHECK TO SEE IF MAXIMUM NUMBER OF ITERATIONS IS EXCEEDED 

IF ( ALPHAT . GT . ANGEND) GO TO 7 
JJJJ=2 

C REENTER CALCULATION AGAIN 
GOTO 6 
7 CONTINUE 

WRITE (6,*) 'MAXIMUM NUMBER OF ITERATIONS EXCEEDED' 
c 

c Uiis section gets and prints the time and date of this run on the outputs 
c 

c CALLGEITIM (IHOUR, IMINUT, ILX3M1,IIXJM2) 

C CALLGETDAT (IYEAR,IMONTH,IDAY) 

C DO 4 IZZZ=2 , 4 

C WRITE (IZZZ, 3 )IMONTH,IDAY,IYEAR, IHOUR, IMINOT 

C 3 FORMAT (' DATE: ' ,12, '/' ,12, '/' ,I4,15X, 'TIME: ',12,':', 12) 

C 4 CONTINUE 

c 

c 

C999 CALL FCLOSE 
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999 CONTINUE 
RETORN 
HID 

SUBROUTINE INITLZ 

OOMMON/XXX/ EQUTVD,RFM,PEXH,PAIR,TAIR, 

1 CR, BORE , STROKE , CONROD, FPISTE , FPISTA, 

2 EVO,EVC,AVO,AVC, 

•3 TVJALL/WHJEL, 

4 ACB, ACF, 

5 ANNA,FMEFM, 

6 NPCWF,DEL1,DEL2,DEL3,AAN(2) , 

7 VALEXH,CDE,NIEXH,WIDTHE,ALPHEX(50) ,FEXH(50) , 

8 VAIATR,CDA,NTAIR,WIDIHA,ALPAIR(50) ,FAIR(50) 

COMMON /GEN /ALFHAT , ANGEND , ANGRES , AIRPO , EXHPO , IMOLS , 

1 APAER,APEXH,AREF,AP(2) , 

1 APN(2) , CYCLE, DALPHA, 

1 DREF , DZ , FREF , 

1 I POWER, IREV,GA,GE, GREF, IUNITL, 

1 IUNITP, IUNITT, IUNITK, IUNITW, IUNITQ, 

1 WAR,PBARAB, 

1 PI , PREF , REVENG , REVREF , RPATR , RPEXH , RP ( 2 ) , 

1 RPN(2) , STEP2 , TEXH , TREF , VREF , WREF , XREF , 

1 Z,ZREF, HAIR, HREF, SCMJLT,NC 

CCWMON/CYI/ACSB,ALHIA,ALFHAE,ANNB,ANNC,AC, 

1 ACN,ACR,ALPHAC, 

1 CALVAL, CARBON, 

1 COEFFA(4) ,COEFFB(4) ,COEFFC(4) ,COEFFD(4) , 

1 COEFFE(4) ,COEFFZ(4) ,CRANK,DCYL,mEAT, 

1 DRC,DVCYL,DWCIN,CWCP(2) ,DWCYL, 

1 FCYL,HEATG, 

1 PURITY, PVEL,PCR,PCYLT, 

1 PORTS, RC,RCN,RCR,SURFC, 

1 TCR,TCYLT, VSWEPT, 

1 VCYL,WCAIR,WCIN, 

1 WCOUT,WCYL,WMW(4) ,V7PCNT(4) ,WCYIR, 

1 WN(4) , WORKGA, WORKGE , XCRA, XSTA, 

1 NCD,WCYLT, CWCINT, 

1 TPIST,THEAD,TEXV, CID 

OCWMON/COffi/AIRFL, EQUTV,TEXHS , SHPS , FZA,HRD, SCASW,TREFF , SCAEFF , 

+ AIRVMX,EXHVMX, PORTED, STOIC, WCO,BLBACK,FUELF,BHP 
******************************************************************** 
OONVERT TO SI UNITS 

******************************************************************** 
11 DCYD=DCYI/3 . 281 
XSTA=XSTA/3 . 281 
XCRA=XCRA/3 . 281 
IF(VALEXH.EQ.O.O) GO TO 15 
DO 14 N=1,NTEXH 

14 FEXH(N)=FEXH(N)/3. 281**2 

15 IF (VAIAER.EQ. 0. 0) GO TO 17 

DO 16 N=1,NTAIR 

16 FAIR (N)=FAIR(N)/3. 281**2 

17 CALLUNITP (IUNITP, 1,PCYIT,PBARAB) 



CALL UNITT ( IUNITT , 1 , TCYLT) 

IF (IUNTIW.EQ.l) GO TO 22 
WFUEL=WFUEI/2 • 205 

22 CALL UNITT ( IUNITT, 1,TWALL) 

CALL UNITT (IUNITT, 1 ,TPIST) 

CALL UNITT ( IUNITT, 1,THEAD) 

CALL UNITT ( IUNITT , 1 , TEXV) 

C ******************************************************************** 
C ENTER ENTHALPY AND VISCOSITY COEFFICIENTS , MOLECULAR WEIGHT 
C JANAF CURVE FIT 400 K TO 2400 K 

C ******************************************************************** 
OOEFFA(1)=3. 00845 
COEFFA(2)=2.94628 
OOEFFA(3)=3.29003 
GOEFFA(4)=3.50953 
OOEFFB ( 1 ) =6 . 1504 IE— 4 
OOEFFB ( 2 ) =1 . 03 6 17E-3 
OOEFFB (3 ) =2.694 17E-3 
OOEFFB (4) =6 . 54488E--4 
COEFFC(l)=— 1.0980 OE— 7 
OOEFFC(2) =— 3 . 39775E— 7 
OOEFFC ( 3 ) =-8 . 66074E— 7 
OOEFFC (4) =9. 4 172 2E-8 
COEFFD ( 1 ) =5 . 68 18 6E-12 
OOEFFD ( 2 ) =4 . 75627E-11 
COEFFD (3) =1 . 11529E— 10 
OOEFFD ( 4 ) =-3 . 4 3 6 6 IE-11 
DO 49 1=1,4 

49 OOEFFE (I)=8.3143*(( COEFFA (I) — 1. ) *TREF +COEFFB ( I ) *TREF**2+ 

+ OOEFFC (I) *TREF**3+COEFFD(I) *TREF**4) 

OOEFFZ (1) =4 . 46E— 7 
COEFFZ (2 ) =5 . 4E— 7 
OOEFFZ (3) =3 . 72E— 7 
OOEFFZ (4) =3. 27E— 7 
WMW(1)=28 . 016 
WMW(2)=32 . 000 
WMW(3)=44 . 010 
WMW(4) =18 . 016 
IMOLS=0 

C SET UP ANGIES FOR EACH CYLINDER 
ALFHAC=EVO— CRANK 

IF (ALFHAC . LT . 0 . ) ALFHAOAIPHAC+180 . * CYCLE 
C SET UP NON-DIMENSIONAL PARAMETERS 
FCYD= (PI/4 . 0) *DCYL**2 

30 IF (VALEXH.EQ. 0. 0) GO TO 31 
DO 32 N=1,NTEXH 

32 FEXH (N) =FEXH (N) /FREF 
GOTO 33 

31 CONTINUE 
ZERO=0.0 

43 CALL STROK (EXHPO, ZERO, EVO,XSTA,XCRA,REVENG,ZREF) 

33 IF (VAIAIR.EQ.0.0) GO TO 34 
DO 35 N=1 , NTAIR 

35 FAIR (N)=FAIR(N) /FREF 
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GOTO 36 
34 ZERO=0.0 

CALL STROK (AIRPO , ZERO, AVO, XSTA, XCRA, REVENG , ZREF) 

36 PVEL=2 . 0*XSTA*XREF*REVREF*REVENG 
WORKGA=0. 

WORKGE=0. 

HEATG=0. 

IMOLS=100 

CALL POWER (1,0) 

37 CONTINUE 
RCR=PCR/PREF 

ACR=SQRT ( 0 . 2 87E3 *GE*TCR) /AREF 

********************************************************************* 

CALCULATE INITIAL CONDITIONS IN CYLINDER 
********************************************************************* 

DST=1. 0 

CALL STROK (XST, DST,EVO, XSTA, XCRA, REVENG, ZREF) 

CALLCYLVOL (FCYL, XST, XSTA, CR,DST,VCYL,DVCYL) 

39 VCYL=XSTA*FCYL 
D7CYL=0.0 
41 CONTINUE 

WCYL=GE*RCR*VCYI/ ( ( (ACR) **2) *GREF) 

WCYLR=WCYL 

ORC=0. 

RC=RCR 
AOACR 
WORKGA=0 . 0 
W0RKGE=0 . 0 
HEATG=0. 0 
WCIN=0 . 0 
SUMIMP=0. 

SUMCAP=0. 

WCOUT=0. 0 
AERNET=0. 

WCAIR=0.0 
DWCYL=0. 0 
RETURN 
END 

SUBROUTINE CYINDR(JJJJ) 

DIMENSION XMOL(4) 

C IMPLICIT REAL ( A-H, O-Z ), INTEGER (I-N) 

OOMMON/XXX/ EQUTVD,REM,PEXH,PAIR,TAIR, 

1 CR , BORE , STROKE , CONROD , FPISTE , FPISTA , 

2 EVO,EVC,AVO,AVC, 

3 TWALL,WFUEL, 

4 ACB,ACF, 

5 ANNA,EMEFM, 

6 NPCWF,DEL1,DEL2,DEL3,AAN(2) , 

7 VALEXH,CDE,NTEXH,WIDIHE,ALFHEX(50) ,FEXH(50) , 

8 VAIAIR,CDA,NTAIR,VTDIUA,ALPAIR(50) , FAIR (50) 

COMMON /GEN /ALFHAT , ANGEND , ANGRES , AIRPO , EXHPO , IMOI5 , 



1 APAIR, APEXH,AREF, AP(2) , 

1 APN(2) , CYCLE, DALPHA, 

1 DREF,DZ,FREF, 

1 IPOWER, IREV,GA,GE,GREF, IUNITL, 

1 IUNITP, IUNITT, IUNITK, IUNITVJ, IUNITQ, 

1 WAR,PBARAB, 

1 PI,FREF,REVENG,REVREF,RPAIR,RPEXH,RP(2) , 

1 RPN(2) ,STEP2,TEXH,TREF,VREF,WREF,XREF, 

1 Z,ZREF, HAIR, HREF, SCMULT,NC 

OCMhK>N/CYI/ACSB, ALFHA, ALFHAE, ANNB, ANNC,AC, 

1 ACN,ACR,ALFHAC, 

1 CALVAL, CARBON, 

1 COEFFA(4) ,COEFFB(4) ,COEFFC(4) ,COEFFD(4) , 

1 COEFFE ( 4 ) , COEFFZ ( 4 ) , CRANK, DCYL, CHEAT, 

1 DRC, DVCYL, CWCIN, CWCP(2) ,EWCYL, 

1 FCYL,HEATG, 

1 PURITY, PVEL,PCR,PCYLT, 

1 PORTS, RC,RCN,RCR,SURFC, 

1 TCR,TCYLT,VSWEPT, 

1 VCYL, WCAIR,WCIN, 

1 WCOUT,WCYL,WMW(4) ,WPCNT(4) ,WCYLR, 

1 WN (4) , WORKGA , WORKGE , XCRA , XSTA , 

1 NCD,WCYLT,DWCINT, 

1 TPIST,THEAD,TEXV, CID 

CCWMON/OCWB/AIRFL, EQUIV, TEXHS , SHPS , FZA , HRD, SCASW,TREFF , SCAEFF , 

+ AIRVMX , EXHVMX , PORTED , STOI C , WCO , BLBACK , FUELF , BHP 
OOMMON/CFURE/WCA , WCHG , DUX 
OOMMON /ENERGY/ UINITL,UFINAL,HOUr 

<XMM0N/TEMP/1VJEV , IWGH , IVJGP , TWGC , TVJCR , TVJPR , TMCP , TWCC , IWCH , TEXHP , 

1 TWEVG , TCI , TC2 , TC3 
OCMMON/INPT2/ 

1 XXEVG , XXEVH , XXPSL, XXPIS , XXPRG , XXHED , XXPCR , XXRSL, XXSLV , 

2 AEXVS , AEXV , AEXST , AHEAD , ASLEEV , APIST1 , APIST2 , APISTR , ARING1 , ARING2 , 

3 CEXV,CHEAD,CSLEEV,CPIST,CRING, 

4 HEXHP,HOILP,HCSL,HCHD,NEXHV,DEXHV, 

5 D1 , D2 , D3 , D4 , D5 , D6 , VCV , APP , ARR , ASS 

OOMMON/HEAT/ NHEAT , QSLEV1 , QPIST1 , QHEAD1 , QEXHV1 , QTOT1 , 

1 QSLEV2 , QPIST2 , QHEAD2 , QEXHV2 , QT0T2 , 

2 HC1 , HC2 , HC3 , HC4 , HC5 , HC6 , HC7 , HC8 , 

3 HCT1 , HCT2 , HCT3 , HCT4 , HCT5 , HCT6 , HCT7 , HCT8 
COMMON/STOP/ ISTOP 

RPS=REVENG*REVREF 
2 OCmTNUE 
ACN=AC 
RCN=RC 
HEX=1 

AFN (NEX) =AP (NEX) 

RFN(NEX) =RP(NEX) 

C SET UP CRANK ANGLE 
101 ALPHAC=ALFHAC+DZ*ZREF*REVENG 
IF(ALFHAC.LT.180.0*CYCLE) GO TO 52 
ALPHAC=ALPHAC-1 8 0 . 0*CYCLE 
52 ALPHA=ALPHAC 
ALFHAE=ALPHA 
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Q ********************************************************************* 

C TEST IF CYLINDER OPEN TO EXHAUST AFTER FIRST REVOLUTION 
C IF SO, START GAS EXCHANGE CALCULATIONS 

C IF NOT, GO TO STATEMENT 12 WHICH INITIALIZES 

C VARIABLES AND RETURNS TO MAIN 

C ********************************************************************* 
CRANK=CRANK— DZ * ZREF*REVENG 
IF (CRANK. IE. 0.0) GO TO 102 

150 DWCP(NEX)=0. 

151 PORTS=0. 0 
GOTO 12 

C SET UP CYLINDER CONDITIONS 
102 RC=RC+DZ *DRC 

WCYL==WCYL+DZ*CWCYL 

Q ********************************************************************* 

C CALCULATE CYLINDER VOLUME , VOLUME CHANGE AND SURFACE AREA 

Q ********************************************************************* 

DVCYLE=0. 0 
DST=1. 0 

CALL STROK (XST,DST, ALPHA, XSTA,XCRA,REVENG,ZREF) 

XHTA=XST 

CALLCYLVOL (FCYL,XST,XSTA,CR,DST,VCYL,DVCYL) 

FPISTE=1 . 0 

CALL SREA ( SURFC , VCYL, FCYL, DCYL, XREF , DREF , 

1 FPISTA, FPISTE) 

C CALCULATE SURFACE AREA OF SIEEVE 
SURFC=SURFC — ( FPISTA+FPISTE) *FCYL 
C CALCULATE SURFACE AREA OF PISTON 
SURFP=FPISTA*FCYL 

C CALCULATE SURFACE AREA OF EXHAUST VALVE (S) 

SURFE=FLOAT (NEXHV) *0 . 7854*DEXHV**2/ ( 39 . 37) **2 
C CALCULATE SURFACE AREA OF HEAD 
SURFH=FPISTE*FCYL-SURFE 
C CALCULATE AC 

AC=SQRT ( GE *RC* VCYL/ (WCYL*GREF) ) 

C ********************************************************************* 

C TEST WHETHER PORTS AND VALVES ARE OPEN, IF SO, SET FLAGS 

C ********************************************************************* 
AIR=10. 

IF(ALFHA.GE.AVC .OR. ALFHA.LE.AVO) AIR=0.0 
EXH=10. 

IF ( ALPHA. GE.EVC .OR. ALPHA. LE.EVO) EXH=0.0 
C ********************************************************************* 
C TEST WHETHER GAS EXCHANGE IS COMPUTE 
C IF SO, FRINT OUT MASS FLOW TOTALS 

C CALCULATE GAS EXCHANGE EFFICIENCIES 

C CALL POWER TO DO COMBUSTION CALCULATIONS 

C ********************************************************************* 
C HAS COMPRESSION STARTED YET? 

IF (ALPHA. LT.ACSB) GOTO 152 
DALFHA=DZ*ZREF*REVENG 

C IS THIS THE FIRST TIME STEP AFTER GAS EXCHANGE? IF SO, DO CALCULATIONS . 

C OTHERWISE, SKIP TO BOTTOM 

IF (ALFHA-ACSB.GT.DALPHA*0.99) GO TO 152 



C IF (ALPHA— ACSB . GE . DALPHA) GO TO 152 

C ********************************************************************* 

C PRINT OUT MASS FLOW TOTALS 

C ********************************************************************* 
R&T ■?=— HFATG— WORKGA+HIN— HOUT— UINITLHJFINAL 
BAL2=BAL2/ (WFUEL*CALVAL) *100. 

WRITE(6, 2223) BAL2 

2223 FOFMAT (2X, 'GAS EXCHANGE ENERGY BALANCE' , F8 . 4 , ' % OF FUEL ENERGY ' ) 

C 

BLBACK=1 . — AIRNET/WCIN 
AIRFL=AIRNET*WREF* 2 . 205*RPS*2 ./CYCLE 
IF (CYCIE . EQ . 4 . ) GO TO 248 
IF(SCMULT.NE. 1. ) GO TO 247 
IF ( PORTED. EQ. 1. ) GO TO 247 

SCASW=WCIN*WREF* . 002871*TAIR/ (PAIR*FCYL*XSTA) 

VEFFYI=1 . -0 . 9811*EXP(— 1 . 1201*SCASW) 

WCAIR=VEFFYI*PAIR*FCYL*XSTA/ (WREF* . 002871*TAIR) 

GOTO 248 

247 WCAIR=WCAIR*SCMULT 
C 

C CALCULATE AND PRINT OUT THE TOTAL AIR SUPPLIED TO THE CYLINDER 
C 

248 WCI=WCIN*WREF*2 . 0*RPS/CYCLE 
CALL UNITY! ( IUNITW , 2 , WCI ) 

IF(IUNITVJ.EQ.2) GO TO 270 
WRITE (4 ,251) WCI 

251 FORMAT (35H TOTAL AIR SUPPLIED PER SECOND(KG)=, Fll. 6) 

GOTO 273 

270 WRITE (4, 252) WCI 

252 FORMAT (35H TOTAL AIR SUPPLIED PER SECOND (LB) =, Fll . 6) 

C 

C CALCULATE AND PRINT OUT THE MASS FICW RATE THROUGH THE EXHAUST VALVE 
C 

273 WOO=4tfOOUT*WREF*2 . 0*RPS/CYCLE 
CALL UNTIW ( IUNITW , 2 , WOO) 

IF(IUNTIW.EQ.2) GO TO 271 
WRITE(4,253)WCO 

253 FORMAT (50H TOTAL MASS FLCW THROUGH EXH VALVE PER SECOND (KG) 

1 Fll. 6) 

GOTO 274 

271 WRITE (4, 254) WOO 

254 FORMAT (50H TOTAL MASS FLCW THROUGH EXH VALVE PER SECOND ( LB) =, 

1 Fll. 6) 

C 

C CALCULATE AND PRINT OUT THE MASS OF AIR RETAINED IN THE CYLINDER 
C 

274 WCA=WCAIR*WREF 

CALL UNTIW (IUNITW, 2, WCA) 

IF ( IUNITW. EQ. 2) GO TO 272 
WRITE (4, 255) WCA 

255 FORMAT ( 34H TOTAL AIR RETAINED PER CYCLE (KG) = , Fll . 6) 

GO TO 275 

272 WRITE (4, 2 56) WCA 

256 FORMAT ( 34H TOTAL AIR RETAINED PER CYCLE (LB) = , Fll . 6) 
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275 SCAVR=WCIN*WREF / ( (PAIR*1. 0E5*VCTL*VREF)/ 

1 (0. 2871*1. 0E3*TAIR)) 

WRITE (4 , 257) SCAVR 

257 FORMAT ( 16H SCAVENGE RATIO= , F6 . 4 ) 

SCASW^WCIN*WREF/ (PAIR*XSTA*XREF*PI* ( (DCYL*EREF) **2)/(4 . 

1*TAIR*. 00287)) 

WRITE(4, 261) SCASW 

261 FORMAT ( 3 3H SCAVENGE RATIO ON SWEPT VOLUME =,F6.4) 

CHEFFY=WCAIR* 100 . 0/WCYL 
RjRITY=CHEFFY/100 . 

one=i.oooo 

FURITY=AMIN1 (RJRITY, ONE) / ( 1 . +WAR) 

CEFFYI=WCAIR*WREF*100 . 0/ ( (PAIR*1 . 0E5*VCYL*VREF) / 

1 (0. 2871*1. 0E3*TAIR) ) 

VEFFYI=CEFFYI*VCYI/ (XSTA*FCYL) 

WRITE (4 , 258) CHEFFY 
WRITE (4 ,259) CEFFYI 
WRITE (4, 260) VEFFYI 

258 FORMAT (2 OH CHARGING EFFICIENCY, 

1 28H BASED ON TRAPPED OONDITION= , F6 . 2 ) 

259 FORMAT (20H CHARGING EFFICIENCY, 

1 43H BASED ON INLET CONDITION AND TRAPPED VOL. = , F6 . 2 ) 

260 FORMAT (2 OH CHARGING EFFICIENCY, 

1 41H BASED ON INLET CONDITION AND SWEPT VOL.=,F6 . 2 ) 

WSUPPLfWCIN * wref 
WCA=WCAIR*WREF 
WCHG=WCYL*WREF 
WCA=AMIN1 (WCA, WCHG) 

WDISP=PAIR*1 . 0E5*VCYL*VREF/ ( . 2871*1 . 0E3*TAIR) 

WDISPS=PAIR*1. 0E5*FCYL*FREF*XSTA*XREF/ ( . 2871*1. 0E3*TAIR) 
TREFF=WCA/WSUPPL* 100 . 

RELCHG=WCHG/WDISP*100 . 

RELCHS=WCHG/WDISPS*100 . 

SCAEFF=WCA/WCHG*100 . 

WRITE ( 4 , 4 0 1 ) TREFF 
WRITE ( 4 , 4 02 ) RELCHG 
WRITE(4,404) RELCHS 
WRITE ( 4 , 4 03 ) SCAEFF 

401 FORMAT (23H TRAPPING EFFICIENCY = , F6 . 2) 

402 FORMAT ( 3 2H RELATIVE CHARGING EFFICIENCY = ,F6.2) 

403 FORMAT ( 2 3H SCAVENGE EFFICIENCY = ,F6.2) 

404 FORMAT (47H RELATIVE CHARGING EFF BASED ON SWEPT VOLUME = ,F6.2) 

WRITE (4, 405) STOIC 

405 FORMAT ( 3 3H STOICHIOMETRIC FUEL AIR RATIO = ,F6.4) 

C ********************************************************************* 
C SET UP PARAMETERS FOR NEXT TIME 

C ********************************************************************* 
WCIN=0 . 0 
WCOUT=0. 0 
SUMLMP=0. 

SUMCAP=0. 

WCAIR=0.0 

AIRNET=0. 
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c ********************************************************************* 
C CALCULATE TRAPPED PRESSURE AND TEMPERATURE 

C ********************************************************************* 
N=1 

PCYLT=FREF* ( RC — ( RC— RCN ) * (ALPHA-ACSB) /DALPHA) 

ACT=AC — ( AC— ACN ) * (ALPHA-ACSB) /DALPHA 
TCYLT= (GREF*TREF*ACT**2 ) /GE 
IF(IPOWER.EQ.O) GO TO 188 
• IM0LS=100 

Q ******** *************************************************** ********** 

C CALL POWER TO DO COMBUSTION CALCULATIONS 

C ********************************************************************* 
CALL POWER (1,1) 

IF(ISTOP.EQ.l) RETURN 
C 

188 RCR=PCR/PREF 

ACR=SQRT (GE*TCR/ (GREF*TREF) ) 

DST=1. 0 

CALL STROK (XST, DST , EVO , XSTA , XCRA, REVENG , ZREF) 

CALL CYLVOL ( FCYL , XST , XSTA , CR , DST , VCYL , DVCYL) 

WCYL==GE*RCR*VCYL/ ( ( (ACR) **2) *GREF) 

WDRKGE=0.0 
W0RKGA=0 . 0 
HEATG=0. 0 
QSLEV2=0. 

QPIST2=0 . 

QHEAD2=0. 

QEXHV2=0. 

H0UT=0. 0 
HIN=0.0 
HOUT2=0 . 0 
HIN2=0.0 
HC5=0. 

HC6=0. 

HC7=0. 

HC8=0. 

HCT5=0. 

HCT6=0. 

HCT7=0. 

HCT8=0. 

C STORE RELEASE MASS 
WCYLR=WCYL 

C TEST IF BOTH AIR AND EXHAUST CLOSED 
152 IF(EXH.NE.O.O) GO TO 124 
IF (AIR. NE. O.O) GO TO 124 
GOTO 150 

C ********************************************************************* 
C CALCULATE GAS EXCHANGE WORK 

C ********************************************************************* 

124 W0RKGA=WORKGA+0 . 5*PREF*VREF* (RC+RCN) * (DVCYL-BVCYIZ) *DZ*100 . 
WORKGE=WORKGE+ 0 . 5*PREF*VREF* (RC+RCN) *DVCYLE*DZ*100 . 

Q ********************************************************************* 

C CALCULATE FIOW THROUGH VALVES/PORTS 

Q ********************************************************************* 
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FOFTS=100.0 

IF(EXH.GT.O.O) GO TO 129 

C EXHAUST VALVE/PORT IS CLOSED, SET MASS FLOW EQUAL TO ZERO 
DWCP(NEX) =0 . 

GOTO 126 
C 

C EXHAUST VALVES/PORTS OPEN 

C CHECK WHETHER ENGINE HAS EXHAUST VALVES OR PORTS 

C 

129 IF (VAIEXH.NE. 0. 0) GO TO 140 

C ENGINE IS PORTED . CALCULATE PORT FLOW AREA AND NONDIMENSIONALIZE . 
XHTE=XHTA 

141 FEXHV=WIDIHE*PI*DCYL* (XHTE-EXHPO) 

FEXHV=FEXHV*CREF*XREF/FREF 
GOTO 542 

C ENGINE IS VALVED, USE SPLINE ROUTINE TO INTERPOLATE USER 

C SUPPLIED EXHAUST VALVE FLOW AREA DATA AT THIS CRANK POSITION 

140 CALL SPLFIT(0, 0,NTEXH, 0, 0. 0, ALFHEX, 0, FEXH, 0,ALFHA,FEXHV, ITRP) 
542 CONTINUE 
533 FEXHV=CDE* FEXHV 

EXHVMX = AMAX1 (EXHVMX, FEXHV) 

C 

C MASS FLOW EXHAUST 
C 

C CHECK WHETHER FLOW IS IN OR OUT OF CYLINDER 

C 

IF(RPEXH.LT.RC) GO TO 590 

CALL MASS FL (FEXHV, RPEXH,RC, APEXH, GREF, GE , CWCP (NEX) ) 

DWCP(NEX)=— EWCP(NEX) 

GOTO 547 

590 CALL MASSFL ( FEXHV, RC,RFEXH, AC, GREF, GE,DWCP (NEX) ) 

547 CONTINUE 
C TEST FOR AIR PORTS 

126 IF(AIR.GT.O.O) GOTO 130 

127 DWCIN=0.0 
GOTO 147 

C 

C INTAKE VALVES/PORTS OPEN 

C CHECK WHETHER ENGINE HAS INTAKE VALVES OR PORTS 
C 

130 IF(VAIAIR.NE.O.O) GO TO 144 

C ENGINE IS PORTED . CALCULATE PORT FLOW AREA AND NONDIMENSIONALIZE . 
FAIRV=WI DTHA*PI * DCYL* (XHTA-AIRPO) 

FAIRV=FAIRV*EREF*XREF/FREF 
GOTO 189 

C ENGINE IS VALVED, USE SPLINE ROUTINE TO INTERPOLATE USER 

C SUPPLIED INTAKE VALVE FLOW AREA DATA AT THIS CRANK POSITION 

144 CALL SPLFn(0,0,NTAIR,0,0.0,ALPAIR,0, FAIR, 0,ALfflA,FAIRV, ITRP) 
189 CONTINUE 
133 FAIRV=CDA*FAIRV 

AIKVMX = AMAX1 (AIRVMX, FAIRV) 

C MASS FLOW AIR 

IF(RPAIR.LT.RC) GO TO 190 

CALL MASSFL (FAIRV, RPAIR,RC,APAIR, GREF, GA,CWCIN) 
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GOTO 147 

190 CALL MASSFL(FAIRV, RC, RPAIR, AC, GREF , GE , EWCIN) 

DWCIN=— OWCIN 
IF(NCD.EQ.l) DWCIN=0. 

C ********************************************************************* 
MASS BALANCE 

********************************************************************* 
147 DWCYD=DWCIN-DWCP(NEX) 

keep track of amount of new air added to cylinder 

ASSUMES CYLINDER CONTENTS ARE PERFECTLY MIXED 

177 IF ( (ALFHA-EVO) — (DZ*ZREF*REVENG) ) 170,171,171 

170 DZT=(ALPHA-EVO)/(ZREF*REVENG) 

GOTO 172 

171 DZT=DZ 

TEST FIOW DIRECTION THROUGH AIR VALVE 

172 IF(DWCIN.GT.O.O) GO TO 174 
OUTFLOW THROUGH AIR VALVE 
TEST FIOW DIRECTION THROUGH EXHAUST VALVE 

IF(DWCP(NEX) .GT.0.0) GOTO 173 
C OUTFLOW AIR VALVE , INFLOW EXH VALVE 
DWEAIR=^7CAIR*EWCIN/WCYL 
WCAIR=WCAIR+DWEAIR*DZT 
GOTO 176 

C OUTFLOW AIR VALVE , OUTFLOW EXH VALVE 

173 DWEAIR=(DWCIN— EWCP(NEX) ) *WCAIR/WCYL 
WCAXR=WCAIR+DWEAIR*DZT 

GOTO 176 

C TEST FIOW DIRECTION THROUGH EXHAUST VALVE 

174 IF (DWCP(NEX) .GT.0.0) GO TO 175 
C INFLOW AIR VALVE , INFLOW EXH VALVE 

DWEAIR=0. 0 

WCAIR=WCAIR+DWCIN*DZT 
WCIN=WCIN+DWCIN*DZT 
GO TO 176 

C INFLOW AIR VALVE, OUTFLOW EXH VALVE 

175 DWEAIR=WCAIR*DWCP(NEX)/WCYL 
WCAIR=WCAIR+ (DWCIN-DWEAIR) *DZT 
WCIN=WCIN+DWCIN*DZT 

C TOTAL MASS FLOW THROUGH EXHAUST VALVE 

176 WOOUT=WCOUTH-EWCP (NEX) *DZT 

C TOTAL MASS FLOWTHROUGH INTAKE VALVE 
AIRNET=AIRNET+DWCIN*DZT 

C ********************************************************************* 
C CALCULATE HEAT TRANSFER 

C ********************************************************************* 
TC= (GREF*TREF*AC**2 ) /GE 
TEXH=TC 

KHO=RC*PREF*l. 0E2/ (0. 287*TC) 

WMOLS=RC*PREF*VCYL*VREF* 1 . 0E2/ (8 . 3143*TC) 

IF(IPOWER.NE.O) GO TO 160 

WN(1)=0.79*WM0LS 

WN(2)=0.21*WM0LS 
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WN(3)=0. 0 
WN(4)=0.0 
GOTO 168 

160 IF(IMOIS.EQ.O) GOTO 164 

C CALCULATE PERCENT BY WT. OF CONSTITUENTS AT END OF POWER CYCLE 
WT=0.0 
DO 165 1=1,4 

165 WT=WIH-WMW (I) *WN (I) 

DO 166 J=1 , 4 

166 WPCNT ( J) =WMW ( J) *WN (J)/WT 
IMOLS=0 

C CALCULATE INSTANTANEOUS NO. OF HOIS . IN GAS EXCHANGE 

164 WN(l)=WREF*(0.767*WCAIR+WPCNr(l) 

1 * (WCYL-WCAIR) ) /WMW ( 1) 

WN ( 2 ) =WREF* ( 0 . 2 3 3 *WCAIR+WPCNT ( 2 ) 

1 * (WCYL-WCAIR)) /WMW (2) 

WN(3)=WREF*WPCNT(3) * (WCYL-WCAIR) /WMW (3) 

WN ( 4 ) =WREF*WPCNT ( 4 ) * (WCYL-WCAIR) /WMW (4) 

C ******************************************************************** 

C CALCULATE HEAT TRANSFER DATA 

C ******************************************************************** 

168 CPT=0 . 0 
DO 161 1=1,4 

161 CPT=CPT+WN ( I ) * (OOEFFA(I)+2 . 0*COEFFB(I) *TC+3.0*COEFFC(I) *TC**2 

1 +4 . 0*C0EFFD(I) *TC**3) 

CPT=8 .3143 *CPT/ (WCYL*WREF) 

SUMMOL=WN ( 1) +WN ( 2 ) +WN (3 ) +WN (4 ) 

RGAS=8 .3143 *SUMMOL/ (WCYL*WREF) 

GE=CPT/ (CPT-RGAS ) 

SUMIMP=SUMIMP+-TC*CPT*CWCP(1) *DZ 

SUMCAP=SUMCAPfCPT*DWCP ( 1 ) *DZ 

TAVG=1 . 8*SUMIMP/SUMCAP 

FRESS=14 . 5*PREF*RC 

TRANK=1.8*TC 

VOL=61032 . 7 *VREF*VCYL 

WRITE (3, 1020) ALPHA, PRESS, TRANK, GE, VOL 

1020 FORMAT (5 (IX, Ell. 5) ) 

XSTIN=XST*39 .372 

C 

C CALCULATE ENTHALPY OF AIR ENTERING CYLINDER 
HTA=0 . 0 
XMOL(l)=.79 
XMOL ( 2 ) = . 2 1 
DO 21 1=1,2 

21 HTA=HTA+XMOL(I) * (8.3143* ( (OOEFFA(I) ) *TAIRfCOEFFB(I) *TAIR**2 
1 -t-COEFPC(I) *TAIR**3+COEFFD(I) *TAIR**4) — OOEFFE(I) ) 

HTA=HTA/28.85 

C CALCULATE ENTHALPY AND INTERNAL ENERGY OF CYLINDER CONTENTS AT 

C THE END OF THE TIME STEP 
HTF=0. 0 
DO 91 1=1,4 

XMOL ( I ) =WN ( I ) /SUMMOL 

91 HTF=HTF+XMOL ( I ) * (8 .3143* ( (COEFFA(I) ) *TC+COEFFB(I) *TC**2 
1 +OOEFFC(I) *TC**3+C0EFFD(I) *TC**4) — €OEFFE(I) ) 
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HTF=HTF*SUMMOI/ (WREF*WCYL) 

BTF=HTF— 6 . 3142*TC*SUMMOL/ (WREF*WCYL) 

C 

EHOUT=DWCP ( 1 ) *HTF*WREF*DZT 
HOUT=HOUT+DHOUT 
EHIN=DWCIN *HTA*WREF* DZT 
HTN=HIN+CHIN 

C CALCULATE COEFFICIENTS FOR GAS VISCOSITY EQUATION 
A=0.0 
B=0.0 

DO 162 1=1,4 

A=WN(I) *COEFFZ (I) *SQPT(WMW(I) )+A 
162 B=«N(I)*SQRT(WMW(I) )+B 

C ********************************************************************* 
C CALCULATE HEAT TRANSFER DURING GAS EXCHANGE 

C ********************************************************************* 
IF(NHEAT.EQ.O) GO TO 777 
IWALLfTVJGC 
TPIST=TWGP 
TKEAD=TWGH 
TEXV=TWEV 

777 VISCTY= ( A*TC* * 0 . 64 5 ) /B 
OONDY=CPT*VISCTY/0 . 7 
REY=RHO*PVEL*DCYL*DREF/VISCrY 
IF (ANNC.EQ. 0) REY=REY*6.18 
H= (ANNA*OONDY*REY**ANNB) / (DCYL*DREF) 

EHEAT=SURFC*H* (TWALL-TC) *VREF*1000 ./AREF 
EHEAT=CHEAT+ (SURFP* (TPIST-TC) +SURFH* (IHEAD-TC) 

1 +SURFE*(TEXV-TC))*H*VREF*1000./AREF 
C TOTAL HEAT TRANSFER DURING GAS EXCHANGE (KJ) 

HEATG=HEATG+EHEAT*DZ * 1 . OE-B 
C 

DQSLEV=SURFC*H* (TWALL-TC) *VREF*1000 ./AREF 
DQPIST=SURFP*H* (TPIST-TC) *VREF*1000 ./AREF 
DQHEAD=SURFH*H* (THEAD-TC) *VREF*1000./APEF 
DQEXHV=SURFE*H* (TEXV-TC) *VREF*1000./AREF 
QSIEV2=QSIEV2+DQSLEV*DZ/1000 . 

QPIST2=QPIST2+DQPIST*DZ/1000 . 

QHEAD2=QHEAD2+DQHEAD*DZ/1000 . 

QEXHV2=QEXHV2+DQEXHV*DZ/1000 . 

QT0T2=QSLEV2+QPIST2-HQHEAD2+QEXHV2 

c 

HC5=HC5+H*SURFC*2 . 

HC6=HC6+H*SURFP*2 . 

HC7=HC7+H*SURFH*2 . 

HC8=HC8+H*SURFE*2 . 

HCT5=HCT5+H*SURFC*TC*2 . 

HCT6=HCT6+H*SURFP*TC*2 . 

HCT7=HCT7+H*SURFH*TC*2 . 

HCT8=HCT8+H*SURFE*TC*2 . 

******************************************************************** 
FIRST LAW OF THERMODYNAMICS 

ENERGY BALANCE AND PRESSURE CHANGE FOR GAS EXCHANGE 
******************************************************************** 
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CVT=€PT-e . 3143*SUMMOI/ (WREF*WCYL) 

DTEMP= (DHEAT/1000 . +WREF* (DWCIN*HTA-EWCP ( 1) *HTF— EWCYL 
+ *ETF) — PC*PREF*DVCYL*VREF*100 . ) / ( WREF*WCYL*CVT ) 
DRC=RC* (ETEMP/TC+DWCYI/WCYL-CVCYI/VCYL) 

12 IF (PORTS. NE. 0.0) GO TO 100 
RC=RCR 
AOACR 
EROO.O 
DWCYL=0.0 
DWCIN=0 . 0 
100 CONTINUE 
1 RETORN 
END 


SUBROUTINE ELIM(AB,N,NP / NDIM) 

DIMENSION AB (NDIM , NP) 

THIS SUBROUTINE SOLVES A SET OF LINEAR EQUATIONS . 

THE GAUSS ELIMINATION METHOD IS USED, WITH PARTIAL PIVOTING. 
MULTIPLE RIGHT HAND SIDES ARE PERMITTED, THEY SHOULD BE SUPPLIED 
AS COLUMNS THAT AUGMENT THE COEFFICIENT MATRIX 
PARAMETERS ARE — 

AB COEFFICIENT MATRIX AUGMENTED WITH R . H . S . VECTORS 
N NUMBER OF EQUATIONS 

NP TOTAL NUMBER OF COLUMNS IN THE AUGMENTED MATRIX 
NDIM FIRST DIMENSION OF MATRIX AB IN THE CALLING PROGRAM. 

THE SOLUTION VECTOR (S) ARE RETURNED IN THE AUGMENTATION 
COLUMNS OF AB. 

REFERENCE: APPLIED NUMERICAL ANALYSIS , 2ND EDITION 
CURTIS F. GERALD P. 132 

BEGIN THE REDUCTION 
NM1=N— I 
DO 35 1=1, NM1 

FIND THE ROW NUMBER OF THE PIVOT ROW. WE WILL THEN 
INTERCHANGE ROWS TO PUT THE PIVOT ELEMENT ON THE DIAGONAL. 

IFVT=I 
IP1=I+1 
DO 10 J=IP1,N 

IF(ABS (AB(IPVT, I) ) .IT. ABS(AB(J,I) ) ) IPVT=J 
10 CONTINUE 

CHECK TO BE SURE THE PIVOT ELEMENT IS NOT TOO SMALL, IF SO 
PRINT A MESSAGE AND RETURN 
IF( (ABS (AB(IPVT,I) ) .IT. l.E-5)) GO TO 99 
NOW INTERCHANGE, EXCEPT IF THE PIVOT ELEMENT IS ALREADY ON 
THE DIAGONAL, DON'T NEED TO. 

IF(IPVT.EQ.I) GO TO 25 
DO 20 JCOL = 1,NP 
SAVE=AB ( I , JCOL) 

AB ( I , JOOL) =AB ( IPVT , JCOL) 

AB ( IPVT , JCOL) =SAVE 
20 CONTINUE 

C NOW REDUCE ALL ELEMENTS BELOW THE DIAGONAL IN THE I-TH ROW. CHECK 



C FIRST TO SEE IF A ZERO ALREADY PRESENT. IF SO, 

C CAN SKIP REDUCTION FOR THAT ROW. 

25 DO 32 JROW=IPl,N 

IF(AB(JROW,I) .EQ.O) GO TO 32 
RATTO=AB(JROW, I)/AB(I, I) 

DO 30 KC0L=IP1,NP 

AB(JROW, KCOL) = AB(JROW,KCOL) — RATIO*AB(I,KOOL) 

30 CONTINUE 
32 CONTINUE 
35 CONTINUE 

C WE STILL NEED TO CHECK A (N,N) FOR SIZE 
IF(ABS (AB(N,N) ) .LT. l.E-5) GO TO 99 
C NOW WE BACK SUBSTITUTE 
NP1=N+1 

DO 50 KC0L=NP1,NP 

AB (N , KOOL) =AB (N , KCOL) /AB (N , N) 

DO 45 J=2,N 
NVBL=NP1-<J 
L=NVBLfl 

VALLJB=AB (NVBL, KCOL) 

DO 40 K=L,N 

VAHJE=VALUE-AB (NVBL, K) *AB (K, KCOL) 

40 CONTINUE 

AB (NVBL, KCOL) =VAIIJE/AB (NVBL, NVBL) 

45 CONTINUE 
50 CONTINUE 
RETURN 

C MESSAGE FOR A NEAR SINGULAR MATRIX 

99 WRITE (6, 100) 

100 FORMAT (2X, 'SOLUTION IS NOT FEASIBLE. A NEAR ZERO PIVOT WAS ENCOUN 
TIERED. ' ) 

RETURN 

END 


SUBROUTINE POWER (N,IRJR) 

DIMENSION VCYLA( 3 60) ,VCYLE(360) ,VCYLP(360) ,PCP(360) ,TCP(360) , 
1 PRANG(360),CANG(2),PCFINr(2),WNI(4),WNF(4) 

C IMPLICIT REAL ( A-H, O— E ), INTEGER (I-N) 

OOMMON/XXX/ EQUIVD, RIM, PEXH, PAIR, TAIR, 

1 CR , BORE , STROKE , CONROD , FPISTE , FPISTA , 

2 EVO,EVC,AVO,AVC, 

3 TWALL,WFUEL, 

4 ACB,ACF, 

5 ANNA,EMEFM, 

6 NPCWF , DELI , DEL2 , DEL3 , AAN ( 2 ) , 

7 VALEXH, CDE,NTEXH, WIETHE, ALFHEX(50) ,FEXH(50) , 

8 VALAIR,CDA,NrAIR,WICTHA,ALPAIR(50) , FAIR (50) 

COMMON/GEN / ALPHAT , ANGEND , ANGRES , AIRPO , EXHPO , IMOLS , 

1 APAIR,APEXH,AREF,AP(2) , 

1 AFN(2) , CYCLE, DALFUA, 
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1 DREF, DZ, FREF, 

1 IPOWER, IREV,GA,GE,GREF, IUNTTL, 

1 IUNITP, IUNTTT, IUNTTK, IUNIW, IUNTTQ, 

1 WAR, PBARAB, 

1 PI,PREF,REVENG,REVREF,RPAXR,RPEXH,RP(2) , 

1 RFN(2) ,STEP2 ,TEXH,TREF,VREF,WREF,XREF, 

1 Z,ZREF, HAIR, HREF, SCMULT,NC 

OCMMON/CYL/ACSB, ALFHA, ALFHAE, ANNB, ANNC, AC, 

1 ACN,ACR,ALPHAC, 

1 CALVAL, CARBON, 

1 COEFFA(4) ,COEFFB(4) ,C0EFFC(4) ,OOEFFD(4) , 

1 C0EFFE(4) ,COEFFZ(4) , CRANK, DCYL, DHEAT, 

1 mC,DVCYL,EWCIN,CWCP(2) ,LWCYL, 

1 ,FCYL,HEATG, 

1 RJRITY,PVEL,PCR,PCYLT, 

1 PORTS , RC, RCN, RCR, SURFC , 

1 TCR,TCYIiT, VSWEPT, 

1 VCYL,WCAIR,WCIN, 

1 WC0UT,WCYL,WMW(4) ,WPCNT(4) ,WCYIR, 

1 WN (4) , WORKGA , WORKGE , XCRA , XST A , 

1 NCD,WCTLT, EWCINT, 

1 TPIST,THEAD,TEXV,CID 

OOMMON/OOMB/AIRFL,EQUIV,TEXHS,SHPS,FZA,HRD,SCASW,TREFF,SCAEFF, 

+ ALRVMX , EXHVMX , PORTED , STOIC , WCO , BLBACK , FUELF , BHP 
<X>MMON/TEMP/TWEV , TWGH , TWGP , TWGC , TWCR , TWPR , TWCP , TWCC , TWCH , TEXHP , 

1 1WEVG , TCI , TC2 , TC3 

OCMMON/HEAT/ NHEAT , QSLEV1 , QPIST1 , QHEAD1 , QEXHV1 , QT0T1 , 

1 QSLEV2 , QPIST2 , QHEAD2 , QEXHV2 , QTOT2 , 

2 HC1 , HC2 , HC3 , HC4 , HC5 , HC6 , HC7 , HC8 , 

3 HCT1,HCT2,HCT3,HCT4,HCT5,HCT6,HCT7,HCT8 
CCWMON/INPT2/ 

1 XXEVG , XXEVH , XXPSL, XXPIS , XXPRG , XXHED , XXPCR , XXRSL, XXSLV , 

2 AEXVS , AEXV , AEXST , AHEAD , ASLEEV , APIST1 , APIST2 , APISTR, ARING1 , ARING2 , 

3 CEXV,CHEAD,CSEEEV,CPIST,CRING, 

4 HEXHP,HOILP,HCSL,HCHD,NEXHV, DEXHV, 

5 D1,D2 ,D3 ,D4 ,D5,D6,VCV ,APP, ARR,ASS 
OOMMON/CV/ QCV,QCV1,QCV2,QCV3 
OCMMON/FRIC/ QFRIC1 , QFRIC2 , QFRIC3 
OCMMON/CFURE/WCA , WCHG , IMIX 
OOMON/ENERGY/ UINITL,UFINAL,HOOT 
COMMON /STOP/ ISTOP 

OOMMON/WHTWY/ PSUM,RSUM, FUELS, TF, DA, X02 
C ******************************************************************** 
C SET UP ANGLES TO START CALCULATION 

C ******************************************************************** 
RPS=REVENG*REVREF 
ANGLA=ACSB 
ANGLE=ACSB 

ANGLP=180 . 0*CYCLE— ACSB+EVO 

ANGLOO.O 

K=1 

PRANG (1)=ACSB 

C SET UP COMBUSTION ANGLES 
ACBB=ACB 
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IF (ACBB. KT. 180.0) ACBB=180. 0*CYCLE+ACBB 
ACFF=ACF 

IF (NPCWF.EQ.O .AND. DEL3.EQ.0.) HRD=15.+115. *EQUIV 
IF (NPCWF.EQ.O) ACFF=ACBB— 18 0 . *CYCLE+HRD 
IF (ACFF.LT. 180.0) ACFF=180. 0*CYCLE+ACFF 
IF (NPCWF.LT. 0.0) HRD=ACFF-ACBB 
C CALCUIATE TRAPPED VOLUME 
DST=0. 0 

CALL STROK (XST , DST , ANGIA , XSTA , XCRA , REVENG , ZREF) 

CALL CYLVOL(FCYL,XST,XSTA,CR,DST,VCYLA(l) ,DVCY) 

FPISTE=1 . 0 

VCYLE(1)=0.0 

VCYLA (1) =VCYIA (1) *VREF 

VCYIE (1) =VCYLE (1) *VREF 

VCYLP(1)=VCYIA(1)+VCYLE(1) 

WORKA=0.0 
W0RKE=0. 0 
HEATRF=0. 0 
DQCR=0.0 
DQCR1=0. 0 
DQCR2=0. 0 
DQCR3=0. 0 
DQCR4=0. 0 
IOOMB=0 
<2iEADl=0. 

QSIEV1=0. 

QPIST1=0. 

QEXHV1=0. 

QCV=0. 

HC1=0 . 

HC2=0 . 

HC3=0 . 

HC4=0 . 

HCT1=0. 

HCT2=0 . 

HCT3-0 . 

HCT4=0. 

C ********************************************************************* 
CALCUIATE TRAPPED GAS COMPOSITION —NO. OF MDIS OF AER(WNA) , 

NO. OF M0L5 OF RESIDUALS (WNR) , INCREASE IN MOLS OF PRODUCTS (WNP) 
********************************************************************* 
VCYLT=VCYLP(1) 

C CALCUIATE TOTAL NUMBER OF MOLES OF GASES IN CYLINDER 
WNTF=PCYLT*VCYLT* 100.0/(8.3143 *TCYLT) 

IF (IFUR.EQ. 0) GO TO 7 
C CALCUIATE MASS OF RESIDUAL 
WRES=WCHG-WCA 

C CALCUIATE TOTAL MASS IN CYLINDER (AIR+RESICUALfFUEL) 

WTOT=WCHG+WFUEL 

C CALCUIATE STOICHIOMETRIC MASS OF AIR REQUIRED TO BURN FUEL 
WABU=32 . * (CARBON/12 . 01+ (1.— CARBON) /4 . 032) *WPUEI/ -233 
ZERO=0.0 

WAIR=WCA+WRES *MAX 1 (WCA-WABU , ZERO) / (V7TOT-WRES ) 

IF (WABU.GT.WCA) WRITE(4 / 8) 
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8 FORMAT (17H *** WARNING *** , 32HBURNED AIR EXCEEDS AVAILABLE AIR) 

7 CONTINUE 

C CALCULATE NUMBER OF MOLES OF AIR IN CYLINDER 
WNA=PURITY*WNTF 

C CALCULATE NUMBER OF MOLES OF RESIDUAL IN CYLINDER 
WNR=WNTF-WNA 

C CALCULATE INCREASE IN NUMBER OF MOLES OF DUE TO COMBUSTION 
C IB FUEL * (MOLES OF PRODUCTS -MDLES OF AIR) /IB FUEL 
WNP=WFUEL* (1 . 0-CARBON) /4 . 032 
C CALCULATE NUMBER OF MDLS OF (1) NITROGEN, (2) OXYGEN, 

C (3) CARBON DIOXIDE, (4) WATER VAPOUR 

WNF ( 1) =0 . 79*WNA* (WNA+WNP+WNR) / (WNA+WNP) 

IF(WNF(1) .LE.O.) WRITE (6,*) / WNA,R3RITY,WNTF,WNF(1) ' ,WNA,FURITY, 

1 WNTF , WNF ( 1 ) 

WNF (4 ) =2 . 0*WNP*WNR/ (WNA+WNP) 

WNF (3 ) =2 . 016*CARBON*WNF (4) / ( 12 . 01* ( 1 . 0— CARBON) ) 

WNF (2 ) =0 . 2 1*WNF ( 1 ) /0 . 79-WNF ( 3 ) -0 . 5*WNF (4 ) 

C TRAPPED MASS IN CYLINDER 
WMF=0. 0 
DO 10 1=1,4 

10 WMF=WMF+WMW(I) *WNF(I) 

TRAF=28 . 85*WNA/WFUEL 
WAIR=WNF(2) *32 ./ . 233 
FURE=^JAIR/WNTF/28 . 85 
WRITE(4,499) PURE 

499 FORMAT ( 14H TRUE FURITY =,F6.4) 

IRMF==WMF/WFUEL 

EAFR=3 2 . *WNF ( 2 ) * 100 . / ( 2 3 . 3 *WFUEL) 

EQUTV=1./ (STOIC*EAFR) 

WRITE(4,4) EQUIV 

4 FORMAT (21H EQUIVALENCE RATIO = ,F6.4) 

IF (EQUIVD.GT.O. ) WFUEL=32 . *WNF(2) *STOIC*EQUIVD/;233 
PSUM=0. 

RSUM=0. 

RJELS=0. 

DUR=0. 

C STORE TRAPPED PRESSURE, TEMPERATURE 

PCF=PCYLT 
PCP(1)=PCF 

CALL UNITP(IUNITP, 2 ,PCP(1) ,PBARAB) 

TF=TCYLT 
TCP(l) =TF 

CALL UNITT ( IUNITT , 2 , TCP ( 1 ) ) 

C ********************************************************************* 
C CALCULATE INTERNAL ENERGY OF TRAPPED CONTENTS IN KJ 

C ********************************************************************* 

ETF=0. 0 
DO 11 1=1,4 

11 ETF=ETF+WNF(I) * (8 . 3143* ( (COEFFA(I) -1. 0) *TF+COEFFB(I) *TF**2 
1 +<X)EFFC(I)*TF**3+OOEFFD(I)*TF**4)-COEFFE(I)) 

UINTTLfETF 

c ********************************************************************* 
C START NEXT STEP 

C ********************************************************************* 



23 K=K+1 

ANGLC=ANGLC+1 . 0 
ANGIA=ANGLA+1 . 0 
ANGLE=ANGLE+1.0 
PRANG (K) =ANGLA 

IF (FRANG(K) .GE. (180.0*CYCLE) ) PRANG (K)=FRANG(K) — 180.0*CYCLE 

ETI=ETF 

PCI=PCF 

TI=TF 

C TEST IF COMBUSTION HAS STARTED 
IF (ANGIA.LE.ACBB) GO TO 12 

C COMBUSTION HAS STARTED —CALCULATE FUEL BURNT IN STEP 
IF (IGOMB.EQ. 0) KO=K— 1 
IC0MB=100 
CANG ( 1) =ANGLA— 1 . 0 
CANG ( 2 ) =ANGLA 
DUR=ANGIA-ACBB 
ONEKL.OOO 
DA=AMIN1 (DUR, ONE) 

C 

C CALCULATE PER CENT FUEL BURNED AT START AND END OF TIME STEP 
C 

DO 25 1=1,2 

C WRITE(6,*) 'CANG(I) , ACBB' ,CANG (I) ,ACBB 
ANG=CANG (I) — ACBB 
IF (ANG) 26,26,27 

26 PCFTNT(I) =0.0 
GOTO 25 

27 IF (NPCWF.LT. 0) GO TO 290 

IF (NPCWF . EQ • 2 ) CALL BURNRT(ANGIA,ACB,ACF,WFUEL,EWFUEL) 

C WRITE (6,*) ANGIA,EWFUEL 
IF (NPCWF. EQ. 2) GO TO 250 
TAU=ANG/HRD 

C ******************************************************************* 
C WIEBE COMBUSTION MODEL 

C ******************************************************************* 
29 IF (NPCWF. EQ. 0) PCFTNT(I)=100. *(1.-EXP(-DEL1*TAU**DEL2) ) 

C ******************************************************************* 
C WATSON COMBUSTION MODEL 

C ******************************************************************* 
BETA=1 . —DELI *EQUIV* *DEL2 / 1 . 5E4**DEL3 
CP1=2.+1 . 25E— 6* ( 1 . 5E4*REVREF*60 . ) **2 . 4 
CP2=5000. 

<3)1=14 . 2/EQUIV** . 644 
CD2=.79*CD1**.25 

IF (NPCWF. GT.O) PCFTNT(I)=100.*(BETA*(1.—(1.-TAU**CP1)**CP2) + 

+ (1.— BETA) * (1.— EXP(— CD1*TAU**CD2) ) ) 

IF (NPCWF. GE.O) GO TO 25 

C ******************************************************************* 

C WHITEHOUSE-WAY COMBUSTION MODEL 

C ******************************************************************* 
290 P02=PCF*WNF(2)/WNTF 
ZERO=0 . 0 

P02=AMAX1 (P02 , ZERO) 
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DWRJELfO. 

C RAMP KIEL INJECTION SCHEDULE ASSUMED 
TAU=DUR/HRD 
AA2=2 . -AAN ( 1 ) *AAN ( 2 ) 

IF(TAU.LE.AAN(2)) FINJ=WFUEL* (AAN ( 1 ) *TAU+ 

+ ( (AA2— AAN(l) )/ (2 . *AAN(2) ) ) *TAU**2) 

IF(TAU.GE. 1. ) GO TO 295 

IF(TAU.GT.AAN(2) ) FINJ=WFUEL* ( (AAN(1)+AA2) *AAN(2)/2.+ 

+ (AA2/ ( 1 . -AAN ( 2 ) ) * (TAU-TAU**2/2 . +AAN ( 2 ) **2/2 . -AAN ( 2 ) ) ) ) 

295 FINJ=AMIN1 ( FINJ , WKJEL) 

IF(DUR.GT.HRD) FINJ=WFUEL 
IF(FINJ.LE.PSUM) GO TO 250 

IREP=DEL1*FINJ** (1.— DEL2) * (FINJ— PSUM) **DEL2*P02**DEL3 

PSUM=PSUM+PREP*DA 

PSUM=AMIN1 ( PSUM , FINJ ) 

REACT=1 . 2E10*P02*EXP ( —1 . 5E4/TF) * (PSUM— RSUM) / (60 . *REVREF*SQRT (TF) ) 
ZERO=0.0 

REACT=AMAX1 (REACT, ZERO) 

RSUM=RSUM+REACT*DA 
IF(RSUM.LE.PSUM) DWFUEL=REACT*DA 
IF (RSUM. GT. PSUM) CWFUEL=PREP*DA 
IF (REACT . EQ . 0 . ) DWFUEL=PSUM-FUELS 
C MODIFIED BURNING RATE MODEL 
GOTO 250 

C ACF1=ACF +360. 

C ACBl=ACB+0 . 8* (ACF1— ACB) 

C IF (ANGIA.LT. ACB) CWFUEL=0. 

C IF ( ANGIA . GT . ACB .AND. ANGIA. LE. ACB1) CWFUEL 
C 1 =WFUEL* ( 2 . / (ACF1— ACB) ) * ( (ANGIA— ACB) / (ACB1-ACB) ) 

C IF ( ANGIA. GT. ACB 1 .AND. ANGIA. LT.ACF1) CWFUEL 
C 1 =WRJEL* ( 2 . / (ACF1— ACB) ) * ( (ACF1— ANGIA) / (ACF1— ACB1) ) 

C IF (ANGIA . GE . ACF1 ) DWFUEL=0 . 

25 CONTINUE 

******************************************************************** 
CALCULATE HEAT RELEASE DUE TO COMBUSTION OF FUEL 
******************************************************************** 
DWFUEL=0. 01* (PCFINT(2) — PCFINT(l) ) *WFUEL 
250 DQF=DWFUEL*CALVAL 
FUELS=FUELS+ CWFUEL 

******************************************************************** 
CALCULATE NEW GAS COMPOSITION 

******************************************************************** 
STORE OLD COMPOSITION 
WNTT=WNTF 
VMI=WMF 
DO 13 1=1,4 
13 WNI(I)=WNF(I) 

C 

X=DWFUEL* ( 1 . 0-CARBON) /4 . 032 
Y=DWFUEL*CARBON/12 . 01 
WNTF=WNTI+X 
WMF=WMI+DWFUEL 
WM=0 . 5* (WMI+WMF) 

WNF(2)=WNI (2) — K— Y 
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WNF(3)=WNI(3)+Y 
WNF(4)=4"JNI (4)+2.0*X 
DO 14 1=1,4 

14 WN(I)=0.5* (WNF (I)+WNI (I) ) 

GOTO 15 

Q ********************************************************************* 

C COMBUSTION HAS NOT STARTED YET 

Q ********************************************************************* 

12 DQF=0 . 0 
DQCR=0 . 0 
WNTI=WNTF 
WMI=WMF 
WM=WMF 
DO 16 1=1,4 
WNI(I)=WNF(I) 

16WN(I)=*JNF(I) 

Q ********************************************************************* 

C CALCULATE CYLINDER VOLUME AND SURFACE AREA 

Q ********************************************************************* 

15 DST=0. 0 

CALL STROK ( XST , DST , ANGLA , XSTA , XCRA , REVENG , ZREF) 

CALL CYLVOL(FCYL,XST,XSTA,CR, DST,VCYIA(K) ,DVCY) 

FPISTE=1 . 0 
VCYLE(K) =0.0 

17 VCYLP(K) =VCYLA (K) +VCYLE (K) 

CALL SREA(SURFC, VCYLP(K) , FCYL,DCYL,XREF,DREF, 

1 FPISTA, FPISTE) 

C CALCULATE SURFACE AREA OF CYLINDER EXPOSED TO COMBUSTION 
SURFC=SURFC— (FPISTA+FPISTE) *FCYL 
SURFC=SURFC*FREF 

C CALCULATE SURFACE AREA OF PISTON (SQUARE METERS) 

SURFP=FPISTA*FCYL*FREF 
C CALCULATE SURFACE AREA OF EXHAUST VALVE (S) 

SURFE=FTOAT (NEXHV) *0 . 7854*DEXHV**2/ (39 . 37) **2 
C CALCULATE SURFACE AREA OF HEAD EXPOSED TO OOMBUSTION CHAMBER 
SURFH=FPISTE*FCYL*FREF — SURFE 
C 

VCYIA (K) =VCYLA (K) *VREF 
VCYLE (K) =VCYLE (K) *VREF 
VCYLP ( K) =VCYLP (K) *VREF 

IF (VCYIP(K) . LT. VCYLP (K—l) ) VOIMLN=VCYLP(K) 

C AVERAGE VOLUME AND DENSITY 

VC=0 . 5* (VCYLP(K) +VCYLP (K-l) ) 

RHO=WIVVC 

C ESTIMATE GAS TEMPERATURE AT END OF STEP 

TF=TI* (VCYLP (K-l) /VCYLP(K) ) **0 . 4+4 . 0*DQF/WM 
C CALCULATE AVERAGE GAS TEMPERATURE 

22 TM=0.5* (TI+TF) 

********************************************************************** 
CALCULATE CONVECTIVE HEAT TRANSFER IN KJ 
USE ANNAND CORRELATION IF ANNC.GT.O 
USE WOSCHNI CORRELATION IF ANNC.EQ. 0 

********************************************************************** 
CALCULATE CP FOR GAS MIXTURE 
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CPT=0.0 
DO 18 1=1,4 

18 CPT=CPT+WN (I ) * (COEFFA ( I) +2 . 0*OOEFFB ( I ) *TMf 3 . 0*CX)EFFC ( I ) *TM**2 

1 +4.0*OOEFFD(I)*TM**3) 

CFT=8 .3143 *CPT/WM 

C CALCULATE VISCOSITY OF GAS MIXTURE 
A=0.0 
B=0.0 

DO 19 1=1,4 
OWN(I)*SQRT(WMW(I)) 

A=A+COEFFZ (I) *C 

19 B=B+€ 

VISCTY= (A*TM**0 . 645) /B 

IF (VTSCTY . IE . 0 . ) WRITE (6,*) 'VISCIY,WN(I)S\VISCTY,WN 
C CA LCU LATE THERMAL CONDUCTIVITY OF GAS MIXTURE (ASSUME PRANDTL =0.7) 
C0NDY=CPT*VISCTY/0 . 7 
REY=RHO*PVEL*DCYL*DREF/VISCrY 

C IF(REY.LE. 0. ) WRITE(6,*) 'REY,RHO,PVEL,DCYL,DREF,VISCrY', 

C 1 REY , RHO , PVEL, DCYL, DREF , VISCTY 
FMOT=PCYLT* ( VCYLT/VCYLP ( K) ) **1 . 32 
TEHM1=2.28*FVEL 
TERM2=0. 

IF ( ANGLA . GT . ACBB) TERM2=3.24E-3*.7854*DCYL**2*XSTA*TCYLT* 

+ (PCF— FM3T) / (PCYLT*VCYLT) 

IF(ANNC.EQ. 0. ) REY=REY * (TERM1+TERM2 ) /PVEL 
H= (ANNA*OONDY*REY**ANNB) / (DCYL* DREF) 

C 

C UNITS: SURFOSQUARE METERS 
C TM = DEG KELVIN 

C RPS = REV PER SECOND 

C DQCY = KT/DEG CA 

HR1=0. 

HR2=0. 

HR3=0. 

HR4=0. 

IF (NHEAT.EQ.O) GO TO 777 

TWALOIWGC 

TPIST=TWGP 

THEAD=TWGH 

TEXV=TWEV 

777 DQCY= (SURFC*H* (TM— TVJALL) )/ (360. 0*RPS) 

DQCY=DQCY+ (SURFP* (TM-TPIST) +SURFH* (TM-THEAD) +SURFE* (TM-TEXV) ) *H/ 

+ (360.*RPS) 

C 

DQCY1=SURFC*H* (TM-TWALL) / (360 . *RPS) 

DQCY2=SURFP*H* (TM-TPIST) / (360. *RPS) 

DQCY3=SURFH*H* (TM-THEAD) / ( 3 60 . *RPS ) 

DQCY4=SURFE*H* (TM-TEXV) /( 360. *RPS) 

IF (IOOMB.EQ.O) GO TO 20 

C ******************************************************************** 
C RADIANT HEAT TRANSFER (AFTER START OF COMBUSTION ONLY) 

C — ANNAND CORRELATION ONLY — 

C ******************************************************************** 
DQCR=SURFC*ANNC* (TM* * 4 -TWALL* *4 ) / ( 3 60 . *RPS ) 
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DQCR=DQCR+ (SURFP* (TM**4-TPIST**4) +SUREH* (TM**4-5HEAD**4 ) 

1 +SURFE* (TM**4— TEXV**4) ) *ANNC/ (360. *RPS) 

C 

HR1=ANNC* (TM**3+TM**2*TWALLfTM*TWALL**2+TWALL**3) 

HR2=ANNC* (TM**3+TM**2*TPIST+TM*TPIST**2+TPIST**3) 

HR3=ANNC* (TM**3+TM**2*THEADfTM*THEAD**2+THEAD**3) 

HR4=ANNC* (TM**3+TM**2*TEXV +TM*TEXV**2 +TEXV**3) 

C 

DQCR1=SURFC*ANNC* (TM**4— 5WALL**4)/ (360.*RPS) 

DQCR2=SURFP*ANNC* (TM**4-TPIST**4)/ (360. *RPS) 

DQCR3=SURPH*ANNC* (TM**4— THEAD**4)/ (360.*RPS) 

DQCR4=SURFE*ANNC* (TM**4-TEXV**4)/ (360.*RPS) 

20 CONTINUE 

Q ********************************************************************* 

C PRESSURE AT END OF STEP 

C ********************************************************************* 
PCF=PCI * VCYLP ( K— 1 ) *TF*WNTF/ (VCYLP (K) *WNTI*TI) 

C ********************************************************************* 
C WORK DONE BY GAS DURING TIME STEP IN KJ 

C ********************************************************************* 
DWORKA=0 . 5* (PCF+PCI) * (VCYIA(K) -VCYIA (K— 1) ) *100 . 0 
DWORKE=0 . 5* (PCF+PCI) * (VCYLE (K) -VCYLE (K-i) ) *100 . 0 
C ********************************************************************* 
C INTERNAL ENERGY AT END OF TIME STEP 

C ********************************************************************* 
ETF=0.0 
DO 21 1=1,4 

21 ETF=ETF+WNF(I)*(8.3143*( (COEFFA(I) -1. 0) *TF+OOEFFB(I) *TF**2 
1 +OOEFFC(I) *TF**3+OOEFFD(I) *TF**4) -COEFFE(I) ) 

UFINAD=ETF 

***** CADCUIATE HEAT TRANSFER FROM CREVICE VOLUME ********** 

CREVICE VOLUME GAS TEMPERATURE — DEG K 
TCV= (ASS*TWGC+APP*TVTCP+ARR*TWCR) / (ASS+APP+ARR) 

C HEAT TRANSFER DUE TO INCOMING HOT GAS — KJ/CRANKANGLE 
DQCV=GE/ (GE— 1 . ) * (TM/TCV-1 . ) *VCV* (PCF-PCI ) *1 . E2 
IF(PCF.LT.PCI) DQCV=0 . 0 

C H-H-++++++++++++++++++-HHH-++-H-+++++++-f H++ ++4++I I I I I t ++-H H-++4-H-+++++ 

C ********************************************************************* 
C FIRST LAW OF THERMODYNAMICS . FOR BALANCE, ERROR=0 

C ********************************************************************* 

c mil t+t I II H H -l I I I I I It I I I M I I H H4 44 -H HH-+++++++++++++++4 +++++++++++ 
ERROR=EnE— ETI-DQF+DQCY+DQCR+EWORKA+EWORKE+DQCV 
DERIV=0.0 
DO 24 1=1,4 

24 DERXV=DERXV+WNF(I) *8.3143* (COEFFA(I) — 1. 0+2.0*COEFFB(I) *TF 
1 +3.0*O0EFFC(I) *TF**2+4 . 0*OOEFFD(I) *TF**3) 

C SET NEW VALUE OF FINAL TEMPERATURE 

ERROR=ERROR/ DERTV 
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TF^TF— ERROR 

C TEST FOR ACCURACY OF ESTIMATE OF TF . IF NOT ACCURATE ENOUGH 
C RETURN TO LABEL 22 WITH BETTER ESTIMATE OF TF AND REPEAT 

C CALCULATION 

IF(ABS (ERROR) .GT. 0.01) GOTO 22 
C 

RGAS=8.3143*(WN(1)+WN(2)+WN(3)4WN(4) )/WM 
GE=CFT/ (CPT-RGAS) 

TEXH=TF 
FRESS=14 . 5*PCF 
TRANK=1.8*TEXH 
VOL=61032.7*VCYLA(K) 

FERCNT=100 . *FUELS/WFUEL 

WRITE(3,1020) PRANG (K) , PRESS, TRANK, GE, VOL 
1020 FORMAT (5 (IX, Ell. 5) ) 

XSTIN=XST*39 . 372 
C ANSWER 

WORKA=WORKA+DWORKA 

WORKE=WORKE+DWORKE 

HEATRF==HEATRF+DQCY+DQCR 

C 

QSIEV1=QSLEV1+DQCY1+DQCR1 

QPIST1=QPIST1+DQCY2+DQCR2 

QHEAD1^HEAD1+DQCY3+DQCR3 

QEXHV1=QEXHV1+DQCY4+DQCR4 

QT0T1=QSLEV1-KJPIST1-K3HEAD1+QEXHV1 

C 

C SUM CREVICE VOLUME HEAT TRANSFER (KJ) 

C 

QCV=QCV+DQCV 

C 

C EVALUATE INTEGRALS USED IN CALCULATING EFFECTIVE GAS TEMPERATURES 
C 

HC1=HC1+ (H+HR1) *SURFC 
HC2=HC2+ (H+HR2 ) *SURFP 
HC3=HC3+ (H+HR3 ) *SURFH 
HC4=HC4+ (H+HR4 ) *SURFE 
HCTl=HCri+ (H+HR1) *SURFC*TM 
HCT2=HCT2+ (H+HR2 ) *SURFP*TM 
HCT3=HCT3+ (H+HR3 ) *SURFH*TM 
HCT4=HCT4+ (H+HR4 ) *SURFE*TM 
C 

TCP(K) =TF 

CALL UNITT(IUNnT,2,TCP(K) ) 

PCP(K)=PCF 

CALL UNITP(IUNITP, 2 ,PCP(K) , PBARAB) 

Q ******************************************************************** 

C TEST IF POWER CYCLE IS COMPLETE. IF NOT REIURN TO IABEL 23 

Q ******************************************************************** 

IF (ANGLO. LT. (ANGLP— 0. 1) ) GO TO 23 

Q ******************************************************************** 

C FINISH POWER CYCLE. STORE RELEASE FRESSURE AND TEMPERATURE 
C AND CALCULATE INDICATED POWER 

Q ******************************************************************** 
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PCR=PCF 

TCR=TF 

WORKP=WORKA+WORKE 

W0RKT=4roRKP+W3RKGA+W0RKGE 

HEATT=HEATRF-HEATG 

CYCLES=RPS*2 . 0/CYCLE 

POWERA=WORKA*CYCLES 

POWERE=WORKE*CYCLES 

FOWEGA=WORKGA*CYCLES 

POWEGE=WORKGE*CYCLES 

POWERT=WORKT*CYCLES 

VSWEPT=XSTA*FCYL*VREF 

HOP= (WORKT/VSWEFT) *0.01 

HEATFLfWRIEL*CALVAL 

CHECK OVERALL ENERGY BALANCE (FOR COMPRESSION, COMB, EXPANSION) 

BALNCE= UFINALr-UINITL^HEATRF+WORKA-HEATFLH-QCV 
PCBALr=BALNCE/HEATFL* 100 . 

WRITE (6, 30) PCBAL 

30 FORMAT (2X, 'COMP-EXP ENERGY BALANCE =',F8. 4,' % OF FUEL ENERGY'/) 
PRINT OUT ANSWERS 

WRITE(4, 116) FRANG(l) ,PCP(1) ,TCP(1) 

116 FORMAT (6H ANGLE, F6. 1,2 3H TRAPPED PRESSURE, F7.1, 

1 26H TRAPPED TEMPERATURE, F7.1) 

WRITE (4, 117) FRANG(KC) ,PCP(KC) ,TCP(KC) 

117 FORMAT (6H ANGLE, F6.1,23H COMPRESSION PRESSURE , F7 . 1 , 

1 26H COMPRESSION TEMPERATURE , F7 . 1) 

1=0 

33 1=1+1 

IF (PCP(I)^CP(I+1)) 33,33,34 

34 WRITE (4, 118) PRANG(I) ,PCP(I) ,TCP(I) 

TMAX=0. 

IMAX=0. 

DO 33 1=1, K 

EMAX=AMAX1 ( PMAX , PCP ( I ) ) 

33 TMAX=AMAX1(TMAX,TCP(I) ) 

WRITE(4, 118) PRANG (I) , PMAX, TMAX 

118 FORMAT (6H ANGLE, F6. 1,2 3H MAXIMUM PRESSURE, F7. 1, 

1 ' MAXIMUM TEMPERATURE', F7.1) 

WRITE(4,119) FRANG(K) ,PCP(K) ,TCP(K) 

119 FORMAT (6H ANGLE, F6. 1,2 3H RELEASE PRESSURE, F7. 1, 

1 26H RELEASE TEMPERATURE, F7.1) 

WRITE(4, 105) TRAF,EAFR, ACBB 

105 FORMAT ( 17H AIR-FUEL RATIO = , F6 . 2 , 

1 3 OH EFFECTIVE AIR-FUEL RATIO =,F6. 2, 

1 24H START OF COMBUSTION = , F7 . 2 ) 

X=VCYLT/VOLMIN 
XX = X 

Y=VCYLP ( K) /VOLMIN 
WRITE(4 , 106) X, Y 

106 FORMAT (3 OH EFFECTIVE COMPRESSION RATIO =,F6.2, 
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1 32H 


EFFECTIVE EXPANSION RATIO = , F6 . 2 ) 
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IF (IUNITK.EQ.2) GO TO 38 
WRITE(4,107) 

107 FORMAT (2 OH INDICATED POWER KW) 

X=POWERA+POWEGA 

WRITE (4, 108) POWERA, POWEGA, X 

108 FORMAT (2 OHAIR POWER CYCIE =,F8.2, 17H GAS EXCHANGE =, F8 . 2 , 

1 10H TOTAL =, F8 . 2) 

39 WRITE (4, 111) POWERT,FMIP 

111 FORMAT ( 16HC0MBINED POWER = , F8 . 2 , 

1 36H MEAN INDICATED PRESSURE (BARS) =,F8. 2) 

GOTO 40 
38 WRITE (4, 112) 

112 FORMAT (23H INDICATED HORSE POWER) 

A=POWERA*l . 34 

B=POWEGA*1.34 

X=A+B 

WRITE(4, 108) A, B,X 
41 A=P0WERT*1.34 
B=PMIP*14 . 5 
AMEP=B 

IF (ALPHAT.LE.540.) GO TO 50 
C 

C CALCULATE EXHAUST TEMPERATURE FROM ENERGY BALANCE 
C 

FUELF=WHJEL*CYCLES*2 . 205 
FZA = FUELF/AIRFL 
YY = HEATT/HEATFL 

DEIH= ( FZA*CALVAL* ( 1 . — YY — WORKT/HEATFL) +HAIR-HREF) / ( 1 . +FZA) 

CALL BAL ( CARBON , FZA , TREF , DEIH , TEXHC , GAM , DUM) 

TEXHC=TEXHC*1 . 8 
C 

YY = YY*100 . 

SHP = A 

IF(ISTOP.NE.l) WRITE(2 , 1000) FUELF,SCASW,TREFF,SCAEFF, TEXHC, 

+ SHP, EURE, AIRFL,EQUIV,FMAX,TMAX,YY 
1000 FOFWAT (F10. 7,T13 , F8. 5,T25,F7. 3 ,T37 , F7 . 3,T47,F7. 1,T56,F7 . 2,T65, 
+F7.4,T77,F8.5,T86,F7.4,T96,F7.1,T107,F7.1,T116,F7.2) 

C 

WRITE(4, 113) A,B 

113 FORMAT (17H COMBINED POWER =,F8 .2 , 

1 38H MEAN INDICATED PRESSURE (P.S.I.) =,F8.2) 

C 

40 A=100 . 0*HEATRF/HEATFL 
B=— 100 . 0*HEATG/HEATFL 
X=100 . 0*HEATT/HEATFL 
WRITE(4,114) A,B,X 

114 FORMAT (35H PERCENT HEAT LOSS POWER CYCIE =,F6.2, 

1 17H GAS EXCHANGE =,F6. 2,1 OH TOTAL=,F6.2) 

C 

Y=100 . 0*W0RKT/HEATFL 
WRITE(4,115) Y 

115 FORMAT ( 3 1H INDICATED THERMAL EFFICIENCY = ,F6.2) 

C 

C 
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CALL HTRANS 


C 

c 

C ******** CONVERGENCE TEST ************************* 

C 

IF (ABS (SHP— SHPS) .GT.SHP* . 001) GOTO 45 
IF (ABS (TEXHC-TEXHS) .GT.TEXHC*.001) GO TO 45 
IF (EQUIVD.GT.O. .AND.ABS(EQUTV— EQUIVD) .GT. .001) GO TO 45 
C 

IF(ISTOP.EQ.l) RETORN 
ISTOP=l 
C 

APIST^. 7854* (3 . 281*DCYL*EREF) **2*144 . 

AIRVMX = AIRVMX*1550 . 15*FREF 
EXHVMX = EXHVMX*1550 . 15*FREF 
EXHAP=EXHVMX/ API ST 
AIRAP=AIRVMX/APIST 
WRITE (2 , 2900) XX,EXHAP, AIRAP,GAM 
2900 FOFMAT (IX, 28HEFFECTIVE COMPRESSION RATIO ,F8 . 4 , 5X, 

+16HEXHVMX/APISTON= , F9 . 4 , 5X , 16HAIRVMX/APIST0N= ,F9.4,5X,4HGE= , 

+ F7.4/) 

WRITE (2, 1010) AMEP , GA , BLRACK , Y 

1010 FOFMAT (IX, 39KINDICATED MEAN EFFECTIVE PRESSURE (PSI) ,F11.4,4X, 

+ 3HGA= , F8 . 5 , 4X , 18HBLOWBACK FRACTION=,F9 . 5, 4X, 

+ 20HINDICATED THERM EFF=,F6.2) 

Q ********************************************************************* 

C ENGINE FRICTION CALCULATION 

Q ********************************************************************* 

VPSTN=FVEL*3 . 281*60. 

FMEP1=12 . 964 
FMEP2=0. 0030476 
FMEP3=0. 65476E— 6 

FMEP=FMEP1+VPSTN* (FMEP2+VPSTN*FMEP3 ) 

FMEP=FMEP*FMEFM 

FHP=FMEP*XSTA*XREF*3 . 281*APIST*REVREF*2 ./(550. *CTCLE) 

BHP=SHP— FHP 
BSFC=3600 . *FUELF/BHP 
C 

WRITE (2 , 1012) VPSTN, FMEP, EHP, BSFC 
1012 F0IMAT(1X, 'PISTON MEAN VEIOCITY (FFM) ' ,F9.2,4X, 

+ 'FRICTION MEAN EFFECTIVE PRESSURE (PSI) ' ,F9.4,4X, 

+ 'BRAKEHP=',F7.2,4X, 'BSFC=' ,F8 .5/) 

45 CONTINUE 
SHPS = SHP 
TEXHS = TEXHC 
50 CONTINUE 
RETURN 
END 
C 

SUBROUTINE MASSFL ( FHOLE , RUP , RDP , AO , GR , G , WFICW) 

C IMPLICIT REAL ( A-H, O-e ) , INTEGER (I-N) 

C SUBROUTINE TO CALCULATE MASS FLOW RATE OF A COMPRESSIBLE 
C FLUID THROUGH A PASSAGE 
C 
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C EHOIE = 

C HUP = RATIO OF UPSTREAM PRESSURE TO REFERENCE PRESSURE 
RDP = RATIO OF DOWNSTREAM PRESSURE TO REFERENCE PRESSURE 
AO = RATIO OF SPEED OF SOUND OF FLUID TO REFERENCE SPEED OF 
SOUND 

GR = REFERENCE VALUE OF GAMMA, THE RATIO OF SPECIFIC HEATS 
G = GAMMA, THE RATIO OF SPECIFIC HEATS 

WFLOW = MASS FLOW RATE OF FLUID THROUGH PASSAGE (NON— DIMENSIONALIZED) 


CALCULATE MACH NUMBER OF FLOW THROUGH PASSAGE TO DETERMINE 
WHETHER FLOW IS CHOKED 

UM=SQRT( (2.0/(G— 1.0) ) * ( ( (RUP/RDP) ** ( (G— 1.0)/G) )— 1.0) ) 
IF(UM.LT.l.O) GOTO 1 

FLOW IS CHOKED 

FM=G*SQRT( (2.0/(G+1.0) )**( (G+1.0)/(G— 1.0) ) ) 

GOTO 2 

FLOW IS NOT CHOKED 

1 FM=G*U1V ( 1 - CH- C (G— 1.0) /2 . 0) *UM**2) ** ( (G+1.0)/ (2 . 0*(G— 1 . 0) ) ) 
CALCULATE MASS FLOW RATE 

2 WFLOW=FHOIE *RUP* FM/ (GR*A0) 

RETURN 

END 

SUBROUTINE HTRANS 

HEAT TRANSFER CALCULATIONS 

(XMMON/XXX/ EQUTVD,RFM,PEXH, PAIR, T AIR, 

1 CR , BORE , STROKE , CONROD , FPISTE , FPISTA , 

2 EVO,EVC,AVO,AVC, 

3 TV7ALL, WRJEL, 

4 ACB,ACF, 

5 ANNA,FMEFM, 

6 NPCWF,DEKL,DEL2,DEL3,AAN(2) , 

7 VALEXH,CDE,NTEXH,WIDLHE,ALPHEX(50) ,FEXH(50) , 

8 VALAER,CDA,NTAIR,WIDIHA,ALPAIR(50) ,FAIR(50) 

COMMON /GEN/AIPHAT , ANGEND , ANGRES , AIRPO , EXHPO , LMOLS , 

1 APAIR,APEXH,AREF,AP(2) , 

1 AFN(2) , CYCLE, DALFHA, 

1 DREF,DZ,FREF, 

1 IPOWER,IREV,GA,GE,GREF,IUNITL, 

1 IUNITP, IUNITT, IUNTTK, IUNITW, IUNITQ, 

1 WAR,PBARAB, 

1 PI,PREF,REVENG,REVREF,RPAIR,RPEXH,RP(2) , 

1 RFN(2) ,STEP2,TEXH,TREF,VREF,WREF,XREF, 

1 Z,ZREF, HAIR, HREF, SCMULT,NC 
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OOMMON/CYI/ACSB , ALPHA , ALFHAE , ANNB , ANNC , AC , 

1 ACN, ACR, ALFHAC, 

1 CALVAL, CARBON, 

1 00EFFA(4) ,C0EFFB(4) ,OOEFFC(4) ,00EFFD(4) , 

1 COEFFE ( 4 ) , COEFFZ ( 4 ) , CRANK, DCYL, CHEAT, 

1 DRC, DVCYL, EWCIN,EWCP(2) ,DWCYL, 

1 FCYL,HEATG, 

1 FURTIY,PVEL,PCR,PCYLT, 

1 PORTS, RC,RCN,RCR,SURFC, 

1 TCR,TCYLT, VSWEPT, 

1 VCYL,WCAIR,WCIN, 

1 WG0UT,WCYL,WMW(4) ,WPCNT(4) ,WCYIR, 

1 WN(4) ,WORKGA,WORKGE,XCRA,XSTA, 

1 NCD,WCYLT,DWCINT, 

1 TPIST,THEAD,TEXV, CID 

OCWMON/CmB/AIRFL, EQUTV, TEXHS , SHPS , FZA, HRD,SCASW,TREFF, SCAEFF, 

+ AIRVMX , EXHVMX , PORTED , STOIC , WCO , BLBACK , FUELF , BHP 
CXM40N/1EMP/TWEV,TWSH,TV?GP,TVJGC,'IWCR,'IVJPR, , IWCP,TWCC,TWCH,TEXHP, 

1 TV7EVG , TCI , TC2 , TC3 

OCWMON/HEAT1/ QF,QFPC,QC,QCPC,A1,A1PC,A2,A2PC,A3,A3PC, 

1 A4,A4PC,QCYL,QCYLPC 
OOMMON/HEAT2/ QEVG , QEVGPC , QEXHP , QEXHPC , 

1 QTC1,QTC1PC,QTC2 / QTC2PC,QTC3 ,QTC3PC,QOUT,QOUTPC 
COMMON /HEAT/ NHEAT , QSLEV1 , QPIST1 , QHEAD1 , QEXHV1 , QT0T1 , 

1 QSLEV2 , QPIST2 , QHEAD2 , QEXHV2 , QT0T2 , 

2 HC1 , HC2 , HC3 , HC4 , HC5 , HC6 , HC7 , HC8 , 

3 HCTl,HCr2,HCr3,HCr4,HCr5,HCr6,HCr7 / HCT8 
OCMMON/INPT2 / 

1 XXEVG,XXEVH,XXPSL,XXPIS,XXPRG / XXHED,XXPCR,XXRSL,XXSLV, 

2 AEXVS , AEXV , AEXST , AHEAD , ASLEEV , APIST1 , APIST2 , APISTR , ARING1 , ARING2 , 

3 CEXV,CHEAD,CSIEEV,CPIST,CRING, 

4 HEXHP,HOILP,HCSL ; HCHD ; 

5 D1 , D2 , D3 , D4 , D5 , D6 , VCV, APP, ARR, ASS 
OCWMON/CV/ QCV, QCV1 , QCV 2 , QCV3 
OCMMON/FRIC/ QFRIC1 , QFRIC2 , QFRIC3 
DIMENSION A (9, 10) ,AFAHR(9) 

RPS=REVENG*REVREF 

C 

C CALCULATE EFFECTIVE GAS TEMPERATURES 
C 

TEG1=HCT1/HC1 
TEG2=HCT2/HC2 
TEG3=HCT3 /HC3 
TEG4=HCT4/HC4 
TEG5=HCT5/HC5 
TEG6=HCT6/HC6 
TEG7=HCT7/HC7 
TEG8=HCT8/HC8 
C 

HIOT1=HC1+HC5 

HP0T2=HC2+HC6 

HT0T3=HC3+HC7 

HK)T4=HC4+HC8 

HITTl=HCTl+HCr5 



HnT2=HCT2+HCT6 
HTTT3=HCr3+HCr7 
HITT4 =HCT4 +HCT8 
C 

TTBGl=HTITl/HrOTl 
TTEG2=HTIT2 /HT0T2 
TTEG3=Hirr3/HTOT3 
TTEG4=HITr4/HrOT4 
C 

C CALCULATE RESISTANCES FROM NETWORK MODEL 
C 

R1=XXEVG/ (AEXVS*CEXV) 

R2=l ./ (HEXHP*AEXV) 

R3=XXEVH/ (AEXST*CEXV) 

R4=XXHED/ (AHEAD* CHEAD) 

R5=l . / (HCHD* AHEAD) 

R6=XXPIS/ (APIST1*CPIST) 

R7=l ./ (HOILP*APIST2 ) 

R8=XXPSI/ (ASLEEV*CPIST) 

R9=XXPCR/ (APISTR*CPIST) 

R10=XXRSL/ (ARING1*CRING) 

R11=XXSLV/ (ASLEEV*CSLEEV) 

R12=XXPRG/ (ARING2*CRING) 

R13=l./ (HCSL*ASLEEV) 

C 

C WRITE(6, *) 'R(S) ',Rl,R2,R3,R4,R5,R6,R7,R8,R9,R10 f Rll,R12,R13 
C 

R14=3 6 0 . /HTOT1 
R1 5=360. /HTOT2 
R1 6=360. /HTOT3 
R17=360 . /HTOT4 
C 

RX1=3 60 . *RPS/HC1 
RX2=360. *RPS/HC2 
RX3=360. *RPS/HC3 
RX4=360 . *RPS/HC4 
RY1=360 . *RPS/HC5 
RY2=360.*RPS/HC6 
RY3=360. *RPS/HC7 
RY4=360.*RPS/HC8 
C 

QX1= (TEG1-TWALL) /RX1 
QX2= (TEG2-TPIST) /RX2 
QX3= (TEG3-THEAD) /RX3 
QX4= (TEG4 -TEXV) /RX4 
QY1= (TEG5-5WALL) /RY1 
QY2= (TEG6-TPIST) /RY2 
QY3= (TEG7-THEAD) /RY3 
QY4= (TEG8-TEXV) /RY4 
C 

C WRITE(6, *) QX1 , QY1 , QSLEV1 , QSLEV2 
C WRITE (6,*) QX2 , QY2 , QPIST1 , QPIST2 
C WRITE(6, *) QX3 , QY3 , QHEAD1 , QHEAD2 
C WRITE (6, *) QX4,QY4,QEXHV1,QEXHV2 
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c 

c 

C HEAT TRANSFER IN KJ/CYCLE 
C 

QCYL1==QSLEV1— QSLEV2 
QCYL2=QPIST1— QPIST2 
QCYL3=QHEAD1— QHEAD2 
QCTL4=QEXHV1-QEXHV2 
C 

QT0T1=QSLEV1+QPIST1+QHEAD1+QEXHV1 

QT0T2=QSLEV24QPIST2+QHEAD2+QEXHV2 

c 

C CHECK CALCULATIONS 
C 

QC1= (TTEGl-TWALL) /R14/RPS 
QC2= ( TTEG2 -TPIST ) /R15/RPS 
QC3= (TTEG3— THEAD) /R16/RPS 
QC4= (TTEG4 -TEW) /R17/RPS 
C WRITE(6,*) 'QCYL1 , QC1 ' , QCYL1 , QC1 
C WRITE(6, *) ' QCYL2 , QC2 ' , QCYL2 , QC2 
C WRITE ( 6 , * ) 'QCYL3 , QC3 ' , QCYL3 ,QC3 

C WRITE(6, *) 'QCYL4 , QC4 ' , QCYL4 , QC4 
QIOTK2CYIA-K5CY^ 

QCIOT=K3Cl'K3C2-l-QC3+QC4 
C WRITE (6,*) 'QTOT,QCTOT' ,QTOT,QCTOT 

qpckjciot/wfuevcalval*ioo . 

WRITE (6, *) 'HEAT LOSS AS A PERCENTAGE OF FUEL INPUT: ' ,QPC 
C 

C **** CALCULATE HEAT GENERATION DUE TO FRICTION **** 

C 

C PISTON/RING FMEP FROM BISHOP'S EQUATION 

C BISHOP, J.N. "EFFECTS OF DESIGN VARIABLES ON FRICTION AND 
C ECONOMY," SAE PAPER 812A, (1964) 

C 

C MEAN PISTON SPEED (M/S) 

STROKE=XSTA 
BORE=DCYL 
CM=2.*STROKE*RPS 
C OIL DENSITY (KG/M**3) 

DENOIL=55 . 31*16 . 018 

C OIL VISOOSITY (KG/M-SEC) AT 200 DEG F 
VISOIL=250.E— 5*1 . 4878 
C REYNOLDS NUMBER 

RE=CM*DENOIL*BORE/VISOIL 
C FMEP, FROM BISHOP (PA) 

FMEPB= (6 . 2E4*CR**0 . 2/RE) *DENOIL*CM**2 
C WRITE (6,*) 'FMEPB,RE,CM' ,FMEPB,RE,CM 
C FRICTION TORQUE (N-M) 

FRICT=FMEPB*BORE*BC®E*STROKE/8 . 

C FRICTION POWER LOSS (KW) 

FRICP=ERICT*2 . *3 . 1415926*RPS/1000 . 

FRICHP=FRI CP/0. 7456 

C WRITE (6,*) 'PISTON/RING FRICTION HP', FRICHP 
C FRICTION HEAT GENERATION TERMS (KW) 
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CFRIC1=0 . 5 
CFRIC2=0 . 5 

QFRIC1=CFRIC1*CFRIC2 *FRICP 
QFRIC2= ( 1 . — CFRIC2 ) *FRICP 
QFRIC3=(1.— CFRIC1) *CFRIC2*FRICP 

***** CREVICE VOLUME CALTUIATIONS ******* 

CREVICE VOLUME HEAT TRANSFER IN KW 

QCV1=QCV*ARR/ (ARR+ASS+APP) *RPS 
QCV2=QCV*ASS/ (ARR+ASS+APP) *RPS 
QCV3=QCV*APP/ (ARR+ASS+APP) *RPS 

GENERATE COEFFICIENT MATRIX FOR SOLUTION OF SYSTEM OF EQUATIONS 

DO 100 1=1,9 
DO 100 J=l,10 
100 A(I,J)=0. 

C 

A(l, 1)=— 1./R13— i./Rll 
A(1,8)=1./R11 
A(l, 10)=— TC2/R13 
C 

A(2,2)=— I./R7— 1./R6 

A(2,6)=l./R6 
A(2, 10)=— TC1/R7 
C 

A(3, 3)=— 1./R5— 1./R4 
A(3,5)=l./R4 
A(3,10)=— TC3/R5 
C 

A(4,4)=— 1./R1— 1./R2— 1./R3— 1./R17 
A(4,5)=l./R3 

A ( 4 , 10 ) =— IVJEVG/R1— TEXHP/R2 -TTEG4/R17 
C 

A(5,5)=— 1./R4— 1./R3— 1./R16 
A(5,3)=l./R4 
A(5,4)=l./R3 
A(5,10)=— 5TEG3/R16 
C 

A(6,6)=— 1./R8— 1./R9— 1./R6— 1./R15 

A(6,8)=l./R8 

A(6,7)=l./R9 

A(6,2)=l./R6 

A(6, 10)=— QFRIC3— TTEG2/R15 
C 

A(7,7)=— 1./R9— 1./R12 
A(7,6)=l./R9 
A(7,9)=l./R12 
A(7, 10)=— QCV3 
C 

A(8,8)=— I./Rll— 1./R10-1./R8— 1./R14 
A(8,1)=1./R11 
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A(8,9)=l./R10 

A(8,6)=l./R8 

A(8, 10)=— QCV2— QFRIC2— 5TEG1/R14 
C 

A(9,9)=— I./R10— 1./R12 
A(9,8)=l./R10 
A(9,7)=l./R12 
A(9 , 10) =— QCV1— QFRIC1 
C 

C SOLVE SYSTEM OF EQUATIONS TO GET METAL TEMPERATURES 
C SOLUTION VECTOR IS LOCATED IN COLUMN 10 OF COEFFICIENT MATRIX 
C TEMPERATURES ARE IN KELVIN 
C 

CALLELIM(A,9,10,9) 

C 

TWCC=A(1, 10) 

TWCP=A(2 , 10) 

TWCH=A(3,10) 

TV7EV=A(4 , 10) 

TWGH=A(5, 10) 

TWGP=A(6, 10) 

TWPfe=A(7 , 10) 

TWGC=A(8 , 10) 

TWCR=A(9, 10) 

C 

DO 200 1=1,9 

200 AFAHR(I)=A(1, 10) *1.8—460. 

C WRITE(6, 1000) (A(I,10) ,AFAHR(I) ,1=1,9) 

1000 FORMAT (2X, 'COOLANT SIDE SLEEVE TEMPERATURE' ,F9.1, ' K',F11.1,' F'/ 

1 2X, 'OIL SIDE PISTON TEMPERATURE' ,5X,F8.1, ' K',F11.1,' F'/ 

2 2X, 'COOIANT SIDE HEAD TEMPERATURE' ,3X,F8.1, ' K',F11.1,' F' / 

3 2X, ' EXHAUST VALVE TEMPERATURE', 7X,F8.1, ' K',F11.1,' F'/ 

4 2X,'GAS SIDE HEAD TEMPERATURE' ,7X,F8. 1, ' K',F11.1,' F'/ 

5 2X, 'GAS SIDE PISTON CROWN TEMPERATURE F7 . 1 , ' K',F11.1, ' F'/ 

6 2X,' PISTON RIM TEMPERATURE', 10X,F8.1,' K',F11.1,' F'/ 

7 2X, 'GAS SIDE SLEEVE TEMPERATURE ',5X,F8.1,' K',F11.1,' F'/ 

8 2X, 'COMPRESSION RING TEMPERATURE ',4X,F8.1,' K',F11.1,' F'/) 

C 

C CHECK ENERGY BALANCE ON HEAT TRANSFER NETWORK 
C 

C — ENERGY IN — 

QF=<2FRIC1+QFRIC2+QFRIC3 
QC=QCV1+QCV2+QCV3 
Al= (TTEG1-TWGC) /R14 
A2= (TTEG2-TWGP) /R15 
A3= (TTEG3-TW3i) /R16 
A4= (TTEG4-TWEV) /R17 
QCYL=QF-KJC+A1+A2+A3+A4 
QFPC=QF/ QCYL* 100 . 

QCPC=QC/QCYL*100 . 

A1FC=A1/QCYL*100 . 

A2PC=A2/QCYL*100. 

A3PC=A3/QCYL*100 . 

A4PC=A4/QCYL*100 . 
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QCYLPC=100. 

C —ENERGY OUT — 

QEVG= (TWEV-WEVG) /R1 

QEXHP= (TWEV— TEXHP) /R2 

CfTCl= (TWCP-TCl) /R7 

QTC2= (TWCC-TC2 ) /R13 

QfTC3= (TWCH-TC3 ) /R 5 

QOUT^EVG4^EXHP4<3TCl-K3frC2+QTC3 

QEVGPC=QEVG/Q0UT*100 . 

QEXHPC=QEXHP/QOUT* 100 . 

QTC1PCK2TC1/QOUT* 100 . 

QTC2PC=<2TC2/Q0UT*100 . 

Qrc3PC==QTC3/Q0UT*100 . 

Q0UTPC=100. 

C PRINT OUT HEAT TRANSFER SUMMARY 
C 

C WRITE (6,*) ' HEAT TRANSFER SUMMARY — KW' 

C WRirE(6,1001) QF , QFPC , QC , QCPC ; 

C 1 Al, A1PC, A2 , A2PC, A3 , A3PC, A4 , A4PC,QCYL,QCYLPC 

1001 FOFMAT (2X,' FRICTION GENERATED HEAT TRANSFER' ,3X,F10. 3, Fll. 2,' %'/ 

1 2X, 'CREVICE VOLUME HEAT TRANSFER 7 , F17. 3, Fll. 2,' %' / 

2 2X,'HEAT TRANSFER TO CYLINDER WALL' ,F15. 3, Fll. 2,' %'/ 

3 2X, 'HEAT TRANSFER TO PISTON' ,4X,F18. 3, Fll. 2,' %'/ 

4 2X, 'HEAT TRANSFER TO HEAD' ,6X,F18.3 ,F11.2, ' %'/ 

5 2X, 'HEAT TRANSFER TO EXHAUST VALVE' ,F15. 3, Fll. 2,' %'/ 

6 2X, 'TOTAL HEAT INPUT TO NETWORK 7 , 3X,F15. 3 ,F11.2, ' %'/) 

C WRITE(6, 1002) QEVG , QEVGPC , QEXHP , QEXHPC , 

C 1 QTC1,QTC1PC,QTC2 ,QTC2PC,QTC3 ,QTC3PC,QOUT,QOUTPC 

1002 FOR4AT(2X, / HEAT TRANSFER TO EXH VALVE GUIDE', 3X,F10. 3, Fll. 2/ %'/ 

1 2X, 'HEAT TRANSFER TO EXHAUST PORT GASES', F10. 3, Fll. 2/ %'/ 

2 2X, 'HEAT TRANSFER TO PISTON COOLING OIL', F10. 3, Fll. 2,' %'/ 

3 2X, / HEAT TRANSFER TO SIEEVE COOLING OIL' ,F10. 3,F11.2, ' %'/ 

4 2X, ' HEAT TRANSFER TO HEAD COOLING OIL ' ,F10.3,F11.2, ' %'/ 

5 2X, 'TOTAL HEAT TRANSFER FROM NETWORK ' ,F10.3,F11.2, ' %') 

C 

RETURN 

END 

C 

SUBROUTINE INHJT2 
OOMMON/INPT2/ 

1 XXEVG,XXEVH,XXPSL,XXPIS,XXPRG,XXHED,XXPCR,XXRSL / XXSLV / 

2 AEXVS,AEXV,AEXST, AHEAD, AS LEEV,APIST1,APIST2,APISTR,ARING1,ARING2, 

3 CEXV,CHEAD,CSLEEV ; CPIST,CRING, 

4 HEXHP,HOILP,HCSL,HCHD,NEXHV, DEXHV, 

5 D1 , D2 , D3 , D4 , D5 , D6 , VCV , APP , ARR , ASS 
OOMMON/XXX/ EQUIVD,RIW,PEXH,PAIR,TAIR, 

1 CR , BORE , STROKE , CONROD , FPISTE , FPISTA , 

2 EVO,EVC,AVO,AVC, 

3 TWALL,WFUEL, 

4 ACB , ACF f 

5 ANNA,FMEFM, 

6 NPCWF , DELI , DEL2 , DEL3 , AAN ( 2 ) , 

7 VAIEXH,CDE,NTEXH,WIDmE,ALPHEX(50),FEXH(50), 

8 VAIATR,CDA,NTAIR < WIDrHA, ALPAIR(50) , FAIR (50) 
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COMMON /GEN /ALFHAT , ANGEND , ANGRES , AIRPO , EXHPO , IMOLS , 

1 APAIR,APEXH,AREF, AP(2) , 

1 AFN(2) , CYCLE, DALFHA, 

1 DREF,DZ,EREF, 

1 I POWER, XREV,GA,GE, GREF, IUNITL, 

1 IUNITP, IUNITT, IUNITK, IUNITW, IUNITQ, 

1 WAR, PBARAB, 

1 PI,PREF,REVENG,REVREF,RPAIR,RPEXH,RP(2) , 

1 RPN(2) ,STEP2,TEXH,TREF,VREF,WREF,XREF, 

1 Z,ZREF, HAIR, HREF, SCMULT,NC 

OCMMON/CYI/ACSB, ALFHA, ALFHAE, ANNB, ANNC, AC, 

1 ACN,ACR,ALPHAC, 

1 CALVAL, CARBON, 

1 COEFFA(4) ,COEFFB(4) ,COEFFC(4) ,COEFFD(4) , 

1 COEFFE ( 4 ) , COEFFZ ( 4 ) , CRANK, DCYL,DHEAT, 

1 DRC,DVCYL,DWCIN,DWCP(2) ,CWCYL, 

1 FCYL,HEATG, 

1 PURITY, PVEL,PCR,PCYLT, 

1 PORTS , RC , RCN , RCR , SURFC , 

1 TCR,TCYLT, VSWEPT, 

1 VCYL, WCAIR, WCIN, 

1 WOOUT,WCYL,WMW(4) ,WPCNT(4) ,WCYLR, 

1 WN (4) ,WORKGA, WORKGE, XCRA,XSTA, 

1 NCD,WCYLT, EWCINT, 

1 TPIST,THEAD,TEXV, CID 

0CWM0N/CO4B/AIRFL, EQUTV , TEXHS , SHPS , FZA , HRD, SCASW,TREFF , SCAEFF , 

+ AIKVMX , EXHVMX , PORTED , STOIC , WCO , BLBACK , FUELF , BHP 
OOMMON/PEMP/IWEV , 1WGH , TWGP , IWGC , TWCR , WPR , IWCP , TWCC , IWCH , TEXHP , 
1 1WEVG , TCI , TC2 , TC3 
PI=3. 141592654 
C BORE AND STROKE IN FEET 
BORE=DCYL 
STROKE=XSTA 

OPEN DATA FILE FOR HEAT TRANSFER INPUT 

OPEN (10, FILE=' DIES ELIO . INP' , STATUS^ OLD 7 ) 

READ IN INITIAL ESTIMATES OF METAL TEMPERA3URES (DEG F) 

READ(10, *) TWEV , TWGH , TWGP, IWGC , IWEVG 
1VJEV=1000. 

1WGH=1000. 

IWGP=1000. 

1WG01000. 

IWEVG=1000. 

READ IN OIL AND COOLANT TEMPERATURES (DEG F) 

READ(10, *) TC1,TC2 ,TC3 
TC1=200. 

TC2=200. 

TC3=200. 
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CONVERT FROM FAHRENHEIT TO KELVIN 

TWEV= (TWEV+4 60 . ) /I . 8 
TW31= (TWGH+4 60 . ) /1 . 8 
TWGP= (TWGP+460 . ) /1 . 8 
TWGC= (TWGC+4 60 . ) /1 . 8 
TVJEVG= (TWEVG+460 . ) /1 . 8 
TCl=(TCl+460.)/1.8 
TC2=(TC2+460. )/l. 8 
TC3=(TC3+460. )/1.8 

READ IN EXHAUST VALVE PARAMETERS AND CALCULATE EXHAUST VALVE AREA 
DEXHV IN INCHES, AEXV IN SQUARE INCHES 
READ(10, *) NEXHV, DEXHV 
NEXHV=4 
EEXHV=1. 

AEXV=FLOAT (NEXHV) *0.785398*DEXHV**2 

READ IN PISTON RING/GROOVE DIMENSIONS (INCHES) 

READ(10,*) D1,D2,D3,D4,D5,D6 

D1=0. 375 

D2=0. 130 

D3=0. 130 

D4=0. 005 

D5=0. 125 

D6=0. 125 

READ IN CONDUCTION PATH LENGTHS (INCHES) 

READ(10, *) XXEVG , XXEVH , XXPSL 
READ(10, *) XXPIS , XXPRG , XXHED 
READ(10, *) XXPCR,XXRSL,XXSLV 
XXEVG=0. 

XXEVH=0. 

XXPSD=0. 

XXPIS=0. 

XXPRG=0. 

XXHED=0. 

XXPCR=0. 

XXRSL=0. 

XXSLV=0 . 2 

XXEVH=0.93 
XXPIS=2 . 00 
XXHED=0.800 
XXSLV=0.750 

CHECK FOR ZEROS, IF SO, USE DEFAULTS 

IF (XXEVG. LE.l.E-6) XXEVG=3. 

IF (XXEVH. LE. l.E— 6) XXEVH=DEXHV/3 . 

IF(XXPSL.I£. l.E— 6) XXPSL=(BORE*12.)/3. 

IF (XXPIS. LE. l.E— 6) XXPIS=0 . 5 
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IF(XXPRG.LE. l.E— 6) XXPRG=D5 
IF(XXHED.LE. l.E— 6) XXHED=0.25 
IF(XXPCR.LE.l.E-6) XXPCR=0.25 
IF(XXRSL.IE. l.E— 6) XXRSL=D6 
IF(XXSLV. IE. l.E— 6) XXSLV=0.25 
C 

C CONVERT FROM INCHES TO METERS 
C 

• XXEVG=XXEVG/39 . 372 
XXEVH=XXEVH/39 . 372 
XXPSD=XXPSI/39 . 372 
XXPIS=XXPIS/39 . 372 
XXPRG=XXPRG/39 . 372 
XXHED=XXHED/39 .372 
XXPCR=XXPCR/39 . 372 
XXRSL=XXRSI/39 . 372 
XXSLV=XXSLV/39 . 372 
C 

C READ IN HEAT TRANSFER AREAS (SQUARE INCHES) 

C 

C READ(10, *) AEXVS, AEXST 
C READ(10, *) ASLEEV , APIST2 
C READ(10, *) APISTR, ARING1,ARING2 
AEXVS=0. 

AEXST=0. 

ASIEEV=0. 

APIST2=0 . 

APISTR=0 . 

ARING1=0 . 

ARING2=0. 

APISTl=FPISTA*PI/4 . *BORE**2 * 144. 
AHEAD=PI/4.*BORE**2*144. — AEXV 
C 

C CHECK FOR ZEROS, IF SO, USE DEFAULTS 
C 

IF (AEXVS . LE . 1 . E— 9) AEXVS=FLOAT (NEXHV) *PI/4 .*(0.25) **2 
IF (AEXST . LE . 1 . E— 9 ) AEXST=FLOAT (NEXHV) *PI*DEXHV*0 . 050 
IF (ASLEEV. LE . 1 . E-9 ) ASLEEV=PI*BORE* (STROKE/ 4 . ) *144 . 

IF (APIST2 . IE . 1 . E— 9 ) APIST2=APIST1 
IF (APISTR. LE.l. E-9) APISTR=PI*BORE*12.*0.5 
IF (ARING1 . LE . 1 . E— 9 ) ARING1=PI*B0RE*12 . *D5 
IF (ARING2 .LE.l. E-9 ) ARING2=PI/4 . * ( ( BORE* 12 . -2 . *D4 ) **2 
1 — (B0RE*12.-2.*D6)**2) 

C 

C CONVERT FROM SQUARE INCHES TO SQUARE METERS 
C 

AEXVS =AEXVS /39. 372**2 
AEXV =AEXV /39. 372**2 
AEXST =AEXST /39 . 372**2 
ASIEEV=ASLEEV/ 3 9 . 372**2 
APISTl=APISTl/39 . 372**2 
APIST2=APIST2/39 . 372**2 
APISTR=APISTF/39 . 372**2 
ARINGl=ARINGl/39 . 372**2 
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ARING2=ARING2/39 . 372**2 
AHEAD = AHEAD/39 . 372**2 

READ IN THERMAL CONDUCTIVITIES — BIU/ (HR-FT-R) 

READ(10, *) CEXV,CHEAD,CSLEEV,CPIST,CRING 
CEXV=24 . 

CHEAD=22 . 

CSLEEV=22. 

CPIST=60. 

CRING=24 . 8 

CONVERT FROM BIU/ (HR-FT-R) TO KW/ (M-K) 

TT=1. 054*3. 281*1. 8/3600. 

CEXV=CEXV*TT 

CHEAD=CHEAD*TT 

CSIEEV=CSLEEV*TT 

CPIST=CPIST*TT 

CRING=CRING*TT 

READ IN CONVECTIVE COEFFICIENTS — BIU/ (HR-FT**2-R) 
READ(10, *) HEXHP,HOILP,HCSL,HCHD 
CONVERT FROM BIU/ (HR-FT**2-R) TO KW/ (M**2-K) 


HEXHP=10. 

H0ILP=900. 

HCSL=200. 

HCHD=200. 

HOILP=900 . 

HCSL=2000. 

HCHD=2000. 

HOILP=900 . 

HCSL=1600. 

HCHD=1600. 

TT=1. 054*3. 281**2*1. 8/3600. 

HEXHP=HEXHP*TT 

HOI LP=HOI LP*TT 

HCSL=HCSL*TT 

HCHD=HCHD*TT 

CONVERT FROM INCHES TO METERS 

Dl=Dl/39 . 372 
D2=D2/39.372 
D3=D3/39.372 
D4=D4/39. 372 
D5=D5/39 . 372 
D6=D6/39 . 372 

CALCULATE CREVICE VOLUME CHARACTERISTICS 
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B0RE=DCYI/3 . 281 

VCV=PI/4 . * (B0RE**2 — (BORE— 2 . *D4) **2) *D1 

1 + PI/4 . * (BORE**2— (BORE— 2 . *D4— 2 . *D2) **2) * (D3— D5) 

2 + PI/4 . * ( (BORE-2 . *D6) * *2 -(BORE-2 . *D4-2 . *D2) **2) *D5 
APP=PI* (BORE-2 . *D4 ) *D1 

1 + PI/4. *( (BORE-2. *D4)**2-(BORE-2.*D4-2.*D2)**2) 

ARR=PI/4 . * (BORE**2 — (BORE— 2 . *D6) **2)+PI* (BORE-2 . *D6) *D5 
ASS=PI*BORE* (D1+D3-D5) 

C 

CLOSE (10, STATUS^ KEEP') 

C 

RETURN 

END 

C 

SUBROUTINE BUKNRT ( ANGIA , ACB , ACF , WFUEL, DWFUEL) 

C MODIFIED BURNING RATE MODEL — SAWTOOTH SHAPE 
ACFl=ACF+360. 

ACBl=ACB+0 . 8* (ACF1-ACB) 

IF (ANGIA.LT. ACB) DWFUEI>0. 

IF ( ANGIA. GT. ACB .AND. ANGIA. LE.ACB1) DWFUEL 
1 =WFUEL* ( 2 . / ( ACF1-ACB) ) * ( (ANGIA— ACB) / ( ACB1-ACB) ) 

IF ( ANGIA. GT.ACB1 .AND. ANGIA. LT.ACF1) DWFUEL 
1 =WFUEL* ( 2 . / (ACF1-ACB) ) * ( (ACF1— ANGIA) / (ACF1-ACB1) ) 

IF ( ANGIA. GE.ACF1) DWFUEL=0. 

RETURN 

END 

SUBROUTINE SPLFIT (JJ , U , NP, OPT , ANGROT , XKX , IXLGR, FY, IYLGR, X , Y , ITRP) 

C 

C SUBROUTINE SPLFIT IS A GENERALIZED SPLINE CURVE FIT PROGRAM 

C 

C THE INPUT PARAMETERS ARE 

C JJ —OVERALL FLOW ITERATION NUMBER 
C U —INTERVAL NUMBER 

C NP —THE NUMBER OF POINTS DESCRIBING THE CURVE 
C (A MINIMUM OF 3 POINTS ARE REQUIRED) 

C OPT —OPTION FOR COORDINATE SYSTEM ROTATION (INTEGER) 

C (USE ONLY FOR REGULAR CARTESIAN COORDINATE SYSTEM) 

C 0 —NO ROTATION 

C 1 —ROTATION (MUST SPECIFY ANGROT) 

C ANGROT —COORDINATE ROTATION ANGLE (DEGREES) 

C (DUMMY VALUE IF OPT=0) 

C XKX (I) —THE NP X— VALUES DEFINING THE CURVE (25 MAX) 

C (MUST BE IN ASCENDING ORDER) 

C IXLGR —LOGARITHMIC X-AXIS INDICATOR 
C 0 —CARTESIAN COORDINATE AXIS 

C 1 —LOGARITHMIC COORDINATE AXIS 

C FY (I) —THE NP Y— VALUES DEFINING THE CURVE 
C (CORRESPONDING TO THE X-VALUES) 

C IYLGR —LOGARITHMIC Y-AXIS INDICATOR 
C O-CARTESIAN COORDINATE AXIS 

C 1 —LOGARITHMIC COORDINATE AXIS 

C NOTE : THE FUNCTION SLOPES AT THE END POINTS ARE EVALUATED FROM THE FIRST 

C TWO AND LAST TWO DATA INRJT POINTS. THESE POINTS MUST THUS BE 


119 



nooooon 


SELECTED TO GIVE A GOOD SLOPE REPRESENTATION 
X —THE X VALUE WHERE THE FUNCTION VALUE IS DESIRED 

THE OUTHJT PARAMETERS ARE 

Y —THE FUNCTION VALUE AT THE DESIRED X-VAHJE 
ITRP —TRAP FDR NON-CONVERGENCE 

REAL MAT (25, 51) ,L(25) ,F'CT(25) ,MKM1,MK,IK,LR(25) 
DIMENSION XKX (25) ,XK(25) ,FY(25) ,F(25) 

DIMENSION XKR(25) ,FR(25) 

INTEGER OPT 

3TRP=0 
N=NP— i 
NP1=N+1 
NP2=N+2 
NP3=N+3 
NM1=N— 1 
NTOT= (NP1+NP2 ) 

DO 5 1=1, NP 
XK(I)=XKX(I) 

F(I)=FY(I) 

5 CONTINUE 

C CHANGE COORDINATES FOR LOGARITHMIC AXES 

IF(IXLGR .EQ. 0)GOTO 210 
DO 200 1=1, NP 

IF(XK(I) .GT. 0.0) GOTO 198 

ITRP=1 

WRITE(6, 606) 

GO TO 1010 
198 CONTINUE 

XK(I)=AU0G10 (XK(I) ) 

200 OONTTNUE 

IF(X .LE. 0. 0)X=1.E— 5 
X=AL0G10 (X) 

210 OONTTNUE 

IF(IYLGR .EQ. OJGOTO 230 

DO 220 1=1, NP 

IF(F(I) .GT. 0. OJGOTO 202 

ITRP=1 

WRITE (6, 608) 

GO TO 1010 
202 OONTTNUE 

F(I)=ALDG10 (F (I ) ) 

220 CONTINUE 
230 CONTINUE 
C 

IF (OPT .EQ. 0) GO TO 15 
ANGROT=ANGROT*3 . 141593/180 . 

CAN=OOS (ANGROT) 

SAN=SIN (ANGRCT) 

C 

DO 10 1=1, NP1 

XKR(I) =XK(I) *CAN + F (I) *SAN 
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FR(I)=F(I) *CAN — XK(I) *SAN 
10 CONTINUE 
15 CONTINUE 
C 

DO 20 1=1, N 

IF (OFT .EQ. 0) GO TO 20 
m(I)=XKR(I+l) -KKR(I) 

20 L(I)=XK(I+1)-KK(I) 

C — - GE T SLOPES AT THE END POINTS 
YPFST= ( F ( 2 ) — F ( 1 ) ) / ( XK ( 2 ) -XK ( 1 ) ) 

YP1ST= (F (NP1) — F (N) )/ (XK(NPl) —XK(N) ) 

IF (OPT .EQ. 1) YPFSTR=(FR(2) — FR(1) )/ (XKR(2) — XKR(l) ) 

IF (OPT .EQ. 1) YPLSTR=(FR(NP1) — FR(N) )/ (XKR(NPl)— XKR(N) ) 

C INITIALIZE THE MATRIX TO ZERO 

DO 30 1=1, NP1 
DO 30 J=1,NP2 
30 MAT(I, J)=0. 

C ADD THE IDENTITY MATRIX TO GET THE OVERALL MATRIX 

DO 40 J=NP3,NTOT 
DO 40 1=1, NP1 
MAT(I, J)=0. 

IF( (J— I) .EQ. (NP1 + 1) )GOTO 35 
GOTO 40 
35 MAT(I, J) =1. 

40 CONTINUE 
C 

MAT(1, 1) =L(l)/3 . 

MAT(l,2)=L(l)/6. 

MAT(1,NP2) = (F(2) — F(l) )/L(l) — ¥PFST 
MAT (NP1 , N) =L (N) /6 . 

MAT (NP1 , NP1) =L(N) /3 . 

MAT (NP1 , NP2 ) =YPI£T — (F (NP1) — F (N) )/L(N) 

IF (OPT .EQ. 1)MAT(1, l)=LR(l)/3. 

IF (OFT .EQ. 1)MAT(1, 2)=LR(l)/6. 

IF (OPT .EQ. 1)MAT(1,NP2) = (FR(2)-FR(1))/IR(1)-YPFSTR 
IF (OPT .EQ. 1) MAT (NP1 , N) =I£(N) / 6 . 

IF(OPT .EQ. l)MAT(NPl,NPl)=LR(N)/3 . 

IF (OPT .EQ. 1)MAT(NP1,NP2)=YPLSTR — (FR(NPl) — FR(N) )/LR(N) 

C 

DO 50 1=2, N 
IM1=I— 1 
IM2=I— 2 
IP1=I+1 

MAT ( I , IM1) =L (IM1 ) /6 . 

MAT(I,I) = (L(IMl)+L(I))/3. 

MAT(I,IPl)=L(I)/6. 

MAT(I,NP2) = (F(IP1) — F(I) )/L(I) -(F(I)-F(IMl) )/L(IMl) 

IF (OFT .EQ. l)MAT(I,IMl)=LR(IMl)/6. 

IF (OPT .EQ. 1)MAT(I,I) = (LR(IM1)+LR(I) )/3. 

IF (OPT .EQ. l)MAT(I,IPl)=LR(I)/6. 

50 IF (OFT .EQ. 1)MAT(I,NP2) = (FR(IP1) -FR(I) )/IE(I) -(FR(I) -FR(IMl) )/IE( 
*IM1) 

C MAKE THE PIVOT ELEMENT THE LARGEST ELEMENT 


NSW=0 



DO 70 J=1,NP1 

IF(J .EQ. NP1)G0 TO 75 

DO 70 I=J,N 

IP=I+1 

IF(ABS(MAT(IP, J) ) .IT. ABS (MAT(J, J) ) )GO TO 70 

NSW=NSWT 

DO 60 JS=l,NTOT 

STOR=MAT(J, JS) 

MAT(J, JS)=MAT(IP, JS) 

MAT(IP,JS)=STOR 
60 OONTINUE 
70 CONTINUE 
75 CONTINUE 

C REDUCE ELEMENTS IN PIVOT COLUMN TO ZERO, EXCEPT PIVOT 

DO 110 J=1 , NP1 
DO 80 IR=1,NP1 

80 FCT(tR)=MAT(IR,J)/MAT(J,J) 

FCT(J)=0. 

DO 100 IZER=1,NP1 
DO 90 JZER=J,NTOT 

MAT ( IZER , JZER) =MAT (IZER, JZER) -FCT ( I ZER) *MAT ( J , JZER) 

90 CONTINUE 
100 CONTINUE 

110 CONTINUE 

C GET THE DETERMINANT 

DET=1.0 
DO 120 K=1 , NP1 
120 DET=DET*MAT (K , K) 

DET=DET* ( (-1 . ) **NSW) 

C TRAP SINGULARITY 

IF (ABS (MAT (NP1 , NP1 ) ) .LT. l.E-7 .AND. ABS(DET) .IT. l.E-7) GOTO 122 
GOTO 125 

122 WRITE ( 6 , 602 ) U , JJ 
ITRP=1 
GO TO 1010 
125 CONTINUE 

C DIVIDE EACH ROW BY IT'S PIVOT 

DO 130 IPIV=1,NP1 
DIV=MAT ( IPIV , IPIV) 

DO 130 JPIV=1 , NTOT 
MAT ( IPIV , JPIV) =MAT ( I PIV , JPIV) /DIV 
130 CONTINUE 

C FOR A GIVEN X, FIND THE XK THAT BRACKET IT AND CALCULATE GENERATED F 

INI>=0 

IF(X .LE. XK(1) )INI>=1 
IF(IND .EQ. 1) Y=F(1) 

IF(IND .EQ. l)GOTO 1000 
IF(X .GE. XK(NPl) )IND=2 
IF(IND .EQ. 2) Y=F(NP1) 

IF (END .EQ. 2) GO TO 1000 
C 

DO 140 1=2, NP1 
IND=0 

IF(X .EQ. XK(I) ) Y=F(I) 
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IF(X .EQ. XK(I))GOTO 1000 

IF(XK(I— 1) .LT. X .AND. XK(I) .GT. X)G0T0 135 
GOTO 140 
135 CONTINUE 
IM1=I— 1 

MKM1=MAT ( IM1 , NP2 ) 

MK=MAT (I , NP2 ) 

XXM=XK(I-i) 

XX=XK(I) 

FK=F(I) 

FKM1=F (I— i) 

LK=XK(I) — XK(I— 1) 

IF (OPT .EQ. 1)XXM=XKR(I— 1) 

IF (OPT .EQ. 1)XX=XKR(I) 

IF (OPT .EQ. 1) FK=FR(I) 

IF (OPT .EQ. 1) FKM1=FR(I— 1) 

IF (OFT .EQ. 1) IK=XKR(I) — XKR(I— 1) 

GOTO 142 
140 CONTINUE 
142 CONTINUE 
C 

IF (OFT .EQ. 0) GO TO 160 

ANGRT1=ANGR0T 

ANGRT2=ANGROT 

ANGKT3=ANGROT 

ANGRT4=ANGROT 

Y1D= (MKM1* (XX-XXM) **3) / (6 . *LK) + (FKM1/IK— LK*MKMl/6 . ) * (XX-XXM) 
Y2 Lf= ( 1 . 0/TAN (ANGRT1) ) *XXM-X/SIN (ANGRT2 ) 

IF(Y1L .GT. Y2L) INDC=1 
IF(Y2L .GT. Y1L) INDC=— 1 
INDCP=INDC 
INDCPI=— INDCP 
DEIXR= (XX-XXM)/ 10. 

XR=XXM 

C 

DO 150 1=1,30 

XR=XR+DELXR 

IF(XR .GT. XX) GO TO 152 

TERM1= (MKM1* (XX-XR) **3 ) / ( 6 . *LK) 

TERM2= (MK* (XR-XXM) **3 ) / ( 6 . *LK) 

TERM3= (FK/IK-MK*W 6 . ) * (XR-XXM) 

TERM4= (FKMl/LK-LK*MKMl/6 . ) * (XX-XR) 

Y1=TERM1+TERM2 +TERM3 +TERM4 

Y2= ( 1 . 0/TAN ( ANGRT3 ) ) *XR-X/SIN (ANGRT4 ) 

IF(Y1 .GT. Y2) INDC=1 
IF(Y2 .GT. Y1)INDC=— 1 
CRIT=ABS (Yl— ¥2) /ABS (Yl) 

IF(CRIT .IE. 1 . E— 4 ) GO TO 155 
IF(INDC .EQ. INDCPI) DELXR=^)ELXR/10 . 

IF(INDC .EQ. INDCPI ) INDCP=INDCPI 
IF(INDC .EQ. INDCPI) INDCPI=— INDCP 
150 CONTINUE 
C 

GOTO 155 
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152 WRITE (6 , 604 ) I , U , JJ 
ITRP=1 
GO TO 1010 
155 CONTINUE 

ANGINV=— ANGROT 
SANI=SIN (ANGINV) 

CANI=COS (ANGINV) 

Y=Y1*CANI — XR*SANI 

C X=XR*CANI + Y1*SANI 

GO TO 1000 
160 CONTINUE 

TERM1= (MKM1* (XX-X) **3 ) / (6 . *LK) 

TERM2= (MK* (X-XXM) **3) / (6 . *LK) 

TERM3= (FK/LK— MK*LK/ 6 . ) * (X-KXM) 

TERM4= (FKM1/IK— LK*MKMl/6 . ) * (XX-X) 

Y=TEPM1+TEPM2+TERM3+TERM4 
C FORMAT STATEMENTS 

602 FORMAT (/ , IX , 4 4HABORT —SINGULAR MATRI X IN S UBROUTINE SPLFIT,/, — 

*1X, 11HIN INTERVAL, 13 , 24H, OVERALL FLOW ITERATION, 13 ,/) 

604 FORMAT (/, IX, 62HABORT —COORDINATE ROTATION CALCULATIONS HAVE NOT C— 
♦ONVERGED IN , 14 , IX , 3 1HITERATIONS IN SUBROUTINE SPLFIT, / , — 

*1X,11HIN INTERVAL, 13, 24H, OVERALL FLOW ITERATION, 13,/) 

606 FORMAT (/, IX, 61HABORT —NEGATIVE NUMBER SPECIFIED IN LOGARITHMIC X 

*AXIS ARRAY,/) 

608 FORMAT (/, IX, 61HABORT —NEGATIVE NUMBER SPECIFIED IN LOGARITHMIC Y 

♦AXIS ARRAY,/) 

1000 CONTINUE 

IF(IXLGR .EQ. 1)X=(10.0)**(X) 

IF(IYLGR .EQ. 1) Y=(10.0) **(Y) 

1010 CONTINUE 
RETURN 
END 
C 

SUBROUTINE UNTTP (IA, IB, PP, PPB) 

C IMPLICIT REAL (A-H,0-£) , INTEGER (I-4J) 

IF(IB.EQ.2) GOTO 1 
IF(IA.EQ.l) GOTO 2 
IF(IA.EQ. 2) GOTO 3 
PP=(PPfPPB)/14.5 
GOTO 2 

3 PP=PP/14 . 5 
2 RETURN 

HF(IA.EQ.l) GO TO 4 
IF(IA.EQ.2) GOTO 5 
PP=(PPfPPB)*14.5 
GOTO 4 
5 PP=PP*14 . 5 

4 RETURN 
END 

C 

SUBROUTINE UNITE (LA, IB, IT) 

C IMPLICIT REAL(A-H,0— 2) , INTEGER ( I -N) 

IF(IB.EQ.2) GOTO 1 
IF(IA.EQ.l) GOTO 2 
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IF(IA.EQ. 2) GOTO 3 
IF(IA.EQ. 3) GOTO 4 
TT=(TT+460.0)/1.8 
GOTO 2 

3 TT=TT+273 . 0 
GOTO 2 

4 TT=TT/1.8 
2 RETURN 

1 IF(IA.EQ.l) GOTO 5 
IF(IA.EQ. 2) GOTO 6 
IF(IA.BQ. 3) GOTO 7 
TT=1.8*TT— 460.0 
GOTO 5 

6 TT=TT— 273 . 0 
GO TO 5 

7 TT=1.8*TT 

5 RETURN 
END 

C 

SUBROUTINE UNTIV(IA, IB, WW) 

C IMPLICIT REAL ( A— H , O-E ), INTEGER (I-N) 

IF(IB.EQ. 2) GOTO 1 
IF(IA.EQ.l) GOTO 2 
WV^WW/2.205 

2 RETORN 

1 IF(IA.EQ.l) GOTO 3 
WW=WW*2. 205 

3 RETORN 
END 

C 

SUBROUTINE STROK(XSTROK, DSTROK ,THETAE, STROKE, CONROD, RPS , ZR) 

C IMPLICIT REAL(A— H,0— E) , INTEGER (I^J) 

OONRAT=2 . 0*CONROD/STROKE 
PI=3. 14159 

THETA=THETAE*PI/ 18 0 . 0 

FNN=SQRT ( CONRAT* * 2 -SIN ( THETA) **2) 

IF(DSTROK.EQ.O.O) GO TO 1 

DSTROK= (PI/360 . 0) *ZR*RPS*STROKE* (SIN (THETA) +SIN (2 . 0*THETA) / 

1 (2.0*FNN) ) 

1 XSTROK=0 . 5*STROKE* ( 1 . 0+OONRAT-OOS (THETA) — FNN) 

RETURN 

END 

C 

SUBROUTINE CYLVOL ( FCYIG , XSTROK , STROKE , CRR, DSTROK, VCYLG, UJCYIG) 
C IMPLICIT REAL (A-H, 0-2 ) , INTEGER (I-N) 

VCYLG=FCYLG* (XSTROK+STROKE/ (CRR— 1 . 0) ) 

D7CYIG=FCYIG*DSTROK 

RETURN 

END 

C 

SUBROUTINE SREA(SAREA,VCYLL,FCYLL,DCYLL,XR,DR,F1,F2) 

C IMPLICIT REAL ( A-H, O-E) , INTEGER (I-N) 

SAREA=4 . 0*VCYLL*XR/ (DCYLL*DR) + (F1+F2 ) *FCYLL 
RETORN 
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END 

C 

SUBROUTINE BAL( CARBON, FZA,T1,DELH,T2, GAM, H3) 

C OVERALL ENERGY BALANCE 
C FZA=OVERALL FUEI/AIR RATIO 
C T=TEMPERATORE, K 
C DELH=ENTHALPY CHANGE, KJ/KG 
DIMENSION X (4) 

C IMPLICIT REAL ( A-fl, O-E ), INTEGER (I-N) 

COMMON/JANAF/ A(4) ,B(4) ,C(4) ,D(4) 

X(l)=. 79 

X(3)=FZA*CARBON*28. 85/12. 01 
X(4)=FZA*(1. -CARBON) *28. 85/2. 016 
X(2)=.21— X(3) — .5*X(4) 

SUtt=X(l)+X(2)+X(3)+X(4) 

X(1)=X(1)/SUM 

X(2)=X(2)/SUM 

X(3)=X(3)/SUM 

X(4)=X(4)/SUM 

IF (DEIH.EQ.O.) GOTO 15 

H1=0. 

DO 5 1=1,4 

5 H1=H1+X(I) *8.3143* (A(I) *T1+B(I) *T1**2+C(I) *T1**3+D(I) *T1**4) 

H2=H1+DELH* (FZA+1 . ) *28 . 85/SUM 
T2=T1+100 . 

15 H2C=0. 

DO 10 1=1,4 

10 H2C=H2C+X ( I ) *8.3143* (A(I) *T2+B(I) *T2**2+C(I) *T2**3+D(I) *T2**4) 
H3=H2C*SU!V(28.85*(1.+FZA) ) 

CP=0. 

DO 20 1=1,4 

20 CP=CPfX(I) *8 .3143* (A(I)+2 . *B(I) *T2+3 . *C(I) *T2**2+4 . *D(I) *T2**3) 
GAM=CP/ (CP-8 . 3143 ) 

IF(DELH.EQ. 0. ) RETORN 
IF(ABS (H2C— H2) .IE.H2*.0001) GO TO 25 
T2=T2+ (H2— H2C) /CP 
GOTO 15 
25 RETORN 
END 
C 

SUBROUTINE POPEN 

CHARACTER* 16 FNAME1 , FNAME2 , FNAME3 , FNAME4 
COMMON /OPEN/ FNAME1 , FNAME2 , FNAME3 , FNAME4 
OPEN (1, FI LE= ' DIESEL1 . INP ' , STATOS= ' OID ' ) 

OPEN (2, FII£=FNAME1 , STATOS= ' NEW ' ) 

OPEN (3, FILE=FNAME2 , STATOS= ' NEW ' ) 

OPEN (4, FIIE=FNAME3 , STATOS= ' NEW ' ) 

RETORN 

END 

C 

SUBROUTINE FCIOSE 
CLOSE (1,STATOS= / KEEP' ) 

CLOSE (2,STATOS= / KEEP / ) 

CLOSE (3,STATOS= / KEEP / ) 
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CLOSE ( 4 , STATUS= ' KEEP ' ) 

RETURN 

END 

C 

SUBROUTINE INFUT 

C IMPLICIT REAL (A-H,0-e) , INTEGER (I-N) 

DIMENSION TITLE (20) 

COMMON/XXX/ EQUIVD,RFM,PEXH,PAIR,TAIR, 

1 CR , BORE , STROKE , CONROD , FPISTE , FPISTA , 

2 EVO,EVC, AVO,AVC, 

3 TWALL,WFUEL, 

4 ACB , ACF , 

5 ANNA, FMEFM, 

6 NPCWF,DEL1,DEL2,DEL3,AAN(2) , 

7 VALEIXH, CDE,NIEXH, WIETHE ,ALFHEX( 50) ,FEXH(50) , 

8 VALAIR,CDA,NTAIR,WIDTHA,ALPAIR(50) ,FAIR(50) 

COMMON /GEN / ALFHAT , ANGEND , ANGRES , ATRPO , EXHPO , IMOL5 , 

1 APAIR, APEXH, AREF, AP(2) , 

1 AFN(2) , CYCLE, DALPHA, 

1 DREF , DZ , FREF , 

1 I POWER, IREV, GA,GE, GREF, IUNITL, 

1 IUNITP, IUNITT, IUNITK, IUNITW, IUNITQ, 

1 WAR, PBARAB, 

1 PI , PREF , REVENG , REVREF , RPAIR , RFEXH , RP ( 2 ) , 

1 RFN(2) , STEP2 , TEXH , TREF , VREF , WREF , XREF , 

1 Z,ZREF, HAIR, HREF, SCMULT,NC 

GOMMON/CYL/ACSB, ALFHA, ALFHAE, ANNB, ANNC, AC, 

1 ACN, ACR, ALPHAC, 

1 CALVAL, CARBON , 

1 OOEFFA(4) ,COEFFB(4) ,COEFFC(4) ,OOEFFD(4) , 

1 COEFFE(4) , COEFFZ ( 4 ) , CRANK, DCYL,DHEAT, 

1 DRC,DVCYL,DWCIN,DWCP(2) ,EWCYL, 

1 FCYL,HEATG, 

1 PURITY, PVEL,PCR,PCYLT, 

1 PORTS, RC,RCN,RCR,SURFC, 

1 TCR,TCYLT, VSWEPT, 

1 VCYL,WCAER,WCIN, 

1 WOOUT,WCYL,WMW(4) ,WPCNT(4) ,WCYIR, 

1 WN(4) ,WORKGA,WORKGE,XCRA,XSTA, 

1 NCD,WCYLT,DWCINT, 

1 TPIST,THEAD,TEXV,CID 

COMMON /OOMB/ AIRFL , EQOIV , TEXHS , SHPS , FZA , HRD , SCASW , TREFF , SCAEFF , 
+ AIRVMX , EXHVMX , PORTED , STOIC , WCO , BLBACK , FUELF , BHP 
COMMON /CHJRE /WCA , WCHG , IMIX 
WRITE(4, 10) 

WRITE (2, 10) 

10 FORMAT ( 'BENSON DIESEL ENGINE MODEL'/) 

C READ(1, 20) TITLE 
C WRITE(4 , 21) TITLE 
C WRITE (2, 21) TITLE 
C 20 FORMAT (20A4) 

C 21 FORMAT (IX, 2 0A4) 

C 


IPCWER=10 



CYCLE=2 . 

ANGEND=7200. 

P0RTED=2. 

OOOLED=2. 

IF (PORTED. EQ. 0. ) PORTED =1. 

C 

IUNITP=2 

IUNTIT=3 

IUNTTL=2 

HJNTIW=2 

IUNITK=2 

IUNITQ=3 

C 

C READ(1, *) EQUIVDjRrW, PEXH, PAIR,TAIR 
C WRITE (2,*) EQUIVD,RFM, PEXH, PAIR, TAIR 
REVENG=RHV60. 

PBARAB=14.7 
FREF=14 . 7 
TREF=518 . 7 
GREF=1. 4 
REVREF=REVENG 
C 

TEXH=2100. 

GE=1. 33 
GA=1 . 4 
WAR=0. 

C 

STEP2=2 . 

C 

C READ(1, *) CR , BORE , STROKE , CONROD , FPISTA 
C WRITE (2,*) CR, BORE, STROKE, CONROD, FPISTA 
DCYLfBC®E/12 . 

XSTA=STROKE/12 . 

XCRA=OCWROD/12 . 

C 

C READ(1,*) EVO , EVC , AVO , AVC 
C WRITE(2 , *) EVO, EVC, AVO, AVC 
ACSB=EVC 

IF (AVC . GT . EVC) ACSB=AVC 
C 

C READ(1,*) IWALL,WFUEL 
C WRITE(2,*) TV3ALL,WFUEL 
FURITY=0.9 
PCYLT=185 . 

TCYLT=1290. 

CRANK=0. 

NCD=0 

C 

C READ(1,*) ACB,ACF 
C WRITE (2,*) ACB,ACF 
TPIST=IWALL 
IHEAD=IWALL 
TEXV=TWALL 
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CALVAL=42515. 

C READ(1, *) ANNA,FMEIM 
C WRITE (2,*) ANNA,FMEFM 
SCMULT=1. 1 
ANNB=0.7 

C 

IF (SCMULT . EQ . 0 . ) SCMULT=1. 

•IF(FMEEM.EQ.O.) FMER*=1. 

C 

CARB0N=0. 85561 

C READ(1,*) NPCWF, DEKL , DEIS , DEL3 , AAN ( 1) , AAN (2) 

C wrtTE (2 * ) NPCWF DELI , DEL2 , DEL3 , AAN ( 1) , AAN ( 2 ) 

C AAN(1K ANN(2) FUEL INJECTION CONSTANTS FOR WHIIEHOUSE-WAY MODEL 

IF (AAN (1) .EQ.O.) AAN(1)=1. 

IF (AAN (2) .EQ.O.) AAN(2)=1. 

IF (DELI .EQ.O.) DELI =6 . 9 

^?^:A(^% 2 'o 1 1 + (1-^^ N )/< 1 - 008M -_ )> * (32 - 0+ 3- 76 * 28 - 

C AVL HEAT RELEASE MODEL USED WHEN NPCWF - 0 , DEL3 0 . 

C WIEBE HEAT RELEASE MODEL USED WHEN NPCWF=0 , DEIS .GT. 0. 

C WHITEHOUSE-WAY MODEL USED WHEN NPCWF .LT. 0 
C WATSON MODEL USED WHEN NPCWF . GT . 0 

IF (NPCWF. GT.O) HRD=125.0 
IF (NPCWF. EQ.O) HRD=DEL3 


C 

C READ(1, *) VAIEXH 
C WRITE (2 , *) VALEXH 

IF(VALEXH.EQ.O.O) GO TO 8 


C 

c 

c 

c 

c 

c 


READ(1, *) CDE,NTEXH 
WRITE (2,*) CDE,NTEXH 


READ(1,*) (ALPHEX(N) ,FEXH(N) ,N=1,NTEXH) 

WRITE(2,191) (ALPHEX(N) ,FEXH(N) ,N-1,NTEXH) 

91 FORMAT(1X,2E10. 3) 

DO 192 N=1,NTEXH 

92 ALFHEX (N) =ALPHEX (N) +EVO 
GOTO 9 


C 

8 CONTINUE 

C READ(1, *) CDE,WIDTHE 
C WRITE (2,*) CDE,WIDIHE 
C 

9 CONTINUE 

C READ(1,*) VALAIR 
C WRITE(2 f *) VALAIR 

IF (VALAIR. EQ. 0.0) GO TO 200 
C 

C READ(1, *) CDA,NTAIR 
C WRITE (2,*) CDA,NTAIR 

C READ(1,*) (ALPAIR(N) ,FAIR(N) ,N=1,NTAIR) 

C WRITE(2 , *) (ALPAIR(N) ,FAIR(N) ,N=1,NTAIR) 


ooo 


DO 193 N=1,NTAIR 
193 ALPATR (N) =ALPATR(N) +AVO 
GO TO 201 

200 CONTINUE 
KEAD(1,*) CDA,WIDTHA 
WRITE (2,*) CDA,WIDIHA 

201 RETURN 
END 

C 

SUBROUTINE HEADER 

CDMMON/XXX/ EQUIVD,RFM,PEXH,PAIR,TAIR, 

1 CR , BORE , STROKE , CONROD , FPISTE , FPISTA , 

2 EVO,EVC,AVO,AVC, 

3 TVJALL,WFUEL, 

4 ACB , ACT , 

5 ANNA, FMEFW, 

6 NPCWF,DEL1,DEL2,DEL3,AAN(2) , 

7 VAIEXH,CDE,NTEXH,WIDTHE,ALFHEX(50) ,FEXH(50) , 

8 VAIAIR, CDAjlTTAIR, WIDTHA,ALPAIR(50) , FAIR (50) 

COMMON /GEN /ALFHAT , ANGEND , ANGRES , AIRPO , EXHPO , IMOIS , 

1 APAIR,APEXH,AREF,AP(2) , 

1 AFN(2) , CYCLE, DALPHA, 

1 DREF,DZ,FREF, 

1 I POWER, IREV, GA,GE, GREF, IUNITL, 

1 IUNITP, IUNITT, IUNITK, IUNIIVJ, IUNITQ, 

1 WAR,PBARAB, 

1 PI, PREF,REVENG,REVREF,RPAIR,RPEXH,RP(2) , 

1 RFN(2) , STEP2 , TEXH , TREF , VREF , WREF , XREF , 

1 Z,ZREF, HAIR, HREF, SCMULT,NC 

OCWMON/CYI/ACSB , ALFHA , ALFHAE , ANNB , ANNC , AC , 

1 ACN,ACR,ALFHAC, 

1 CALVAL, CARBON, 

1 COEFFA(4) ,OOEFFB(4) ,OOEFFC(4) ,COEFFD(4) , 

1 COEFFE ( 4 ) , COEFFZ ( 4 ) , CRANK, DCYL,DHEAT, 

1 DRC,DVCYL,DWCIN,DWCP(2) ,CWCYL, 

1 FCYL,HEATG, 

1 FURITY,FVEL,PCR,PCYLT, 

1 PORTS, RC,RCN,RCR,SURFC, 

1 TCR,TCYLT,VSWEPT, 

1 VCYL,WCAIR,WCIN, 

1 WCXXJT,WCYL,WMW(4) ,WPCNT(4) ,WCYIR, 

1 WN(4) , WORKGA , WORKGE , XCRA , XSTA , 

1 NCD,WCYLT,DWCINT, 

1 TPIST,THEAD,TEXV,CID 

OOMMON/OOMB/AIRFL, EQUIV, TEXHS , SHPS , FZA, HRD, SCASW,TREFF, SCAEFF, 
+ AIRVMX, EXHVMX, PORTED, STOIC, WOO, BUBACK, FUELF, HIP 
OCWMON/CPURE/WCA , WCHG , IMEX 
WRITE (3, 2048) 

2048 FORMAT (4X, 'ANGLE 7 ,5X, 'PRESSURE' ,3X, # TEMPERATURE 7 ,4X, 

1 ' GAMMA', 5X, 'VOLUME' ,5X, 'MEAN TEMP'/) 

C 

NC=1 
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WRITE (2, 1027) 

WRITE (2, 1028) TWALL,WFUEL, ACB, ACF,TPIST,THEAD,WAR 

1027 FORMAT (43X, 'WALL' ,7X, 'FUEL' ,8X, 'COMBUSTION' , 5X, 'PISTON' , 4X, 

+ 'CHAMB'/ 11X, 13X, 13X, 6X, 'TEMP' ,5X, 

+ 'PER CYCLE ',4X, 'START FINISH', 5X, 'TEMP TEMP' ,5X, 'WAR'/) 

1028 FORMAT (39X, F9.1,E13.3,F9.1,F7.1,F10.1,F9.1,3X,F6.4) 

1035 FORMAT (14H CYLINDER BORE, F8 . 4/7H STROKE, F9. 4/ 

1 50H SURFACE AREA PISTON / CYLINDER CROSS SECTION AREA, F8 . 4 

2 /22H CONNECTING ROD LENGTH, F9. 4) 

DCYL12=DCYL* 12 

XSTA12=XSTA* 12 

XCRA12=XCRA* 12 

CII>=. 7854*DCYL12**2*XSTA12 

WRITE (2, 1035) DCYL12 , XSTA12 , FPISTA, XCRA12 

WRITE(2,1038) CR,CID 

1038 FORMAT ( 2 6H NOMINAL CONCESSION RATIO, F8 .4, 5X, 'SWEPT VOLUME (IN3) ' 
+,F8.3) 

WRITE(2, 1040) EVO,EVC,AVO,AVC,ACSB,AAN(l) ,AAN(2) 

1040 FORMAT (/12H PORT TIMING//4H EVO , F7 . 1/4H EVC, F7 . 1/4H AVO , F7 . 1/ 

1 4H AVC , F7 . 1/32H ANGIE COMPRESSION STROKE BEGINS , F7 . 1 , 5X , 

l'AAN(l)= ' , F6 . 3 , 5X , ' AAN ( 2 ) = ',F6.3) 

IF (VALEXH.GT.5.) GO TO 85 
WRITE (2, 1042) WIDTHE 

1042 FORMAT (28H PERCENTAGE PORT WIDTH (EXH) ,F9.4) 

GOTO 86 

85 WRITE (2, 1043) CDE 

1043 FORMAT (26HCOEFFICIENT CD (EXH) VALVE, F7. 4) 

86 IF (VALAIR. GT. 5. ) GO TO 87 
WRITE(2, 1044) WIDTHA 

1044 FORMAT ( 2 7HPERCENTAGE PORT WIDTH (AIR) ,F9.4) 

GOTO 88 

87 WRITE (2, 1045) CDA 

1045 FORMAT (25HOOEFFICIENT CD (AIR) PORT,F7.4) 

88 SPEED=REVENG*60. 

WRITE (2, 2047) SPEED, PAIR, TAIR,PEXH 
2047 FORMAT (11HSPEED (RFM) ,F10.2,5X,26HSUPPLY AIR PRESSURE (PSIA) , 

+ F8 . 2 , 5X, 19HSUPPLY AIR TEMP (R) , F8 . 2 , 5X , 2 3HEXHAUST PRESSURE (PSIA) 

+ , F8.2/) 

WRITE(2, 1049) CALVAL, CARBON , DELI , DEL2 , DEL3 
1049 FORMAT (22HCAL0RIFIC VALUE (KJ/KG) , F10 . 2 , 4X, 

1 15HCARBON FRACTION, F8. 5, 4X,5HDEL1 ,F6.4,4X,5HDEL2 ,F6.4, 

+4X, 5HDEL3 ,F6.4) 

WRITE (2, 1051) ANNA,IMEFM,ANNB 
1051 FORMAT (1H ,8HANNAND A,F7.4,49X, 

+ 5X, 'FMEFM= ' ,F5.3,/9H ANNAND B,F7. 4) 

WRITE (2, 1052) ANNC 
WRITE (2, 3000) 

3000 FORMAT (/IX, 9HRJEL FLOW, 2X, 9HSCAV RAT. , 3X, 8HTRAPPING , 3X , 

+ 8HSCAVENGE, 2X, 8HTEMP OUT, 3X, 5HP0WER,5X, 4HTRUE,5X, 8HAIR FLOW, 

+ 3X, 5HEQUTV, 5X, 5HPMAX,6X, 5HTMAX,4X, 7HPERCENT) 

WRITE (2, 3010) 

3010 FORMAT (2X, 6HIB/SEC,5X, 6HSWVOL,7X, 3HEFF,8X, 3HEFF,9X, 1HR,8X, 

+ 2HHP, 5X, 6HFURITY, 5X, 6HLB/SEC, 15X, 4HPSIA,8X, 1HR,5X, 

+ 9HHEAT DOSS/) 



oooooo 


1052 F0FMAT(1H , 8HANNAND C,E12.4,2X,16H KW/ (M**2) (K**4)) 

RETURN 

END 

C 

BLOCK DATA 

FOURTH ORDER POLYNOMIAL CURVE FIT TO JANAF TABLE DATA 
IN THE RANGE 400 K TO 2400 K 

THIS DATA CAN BE USED TO FIND ENTHALPY, INTERNAL ENERGY , AND 
SPECIFIC HEAT USING THE FOLLOWING EQUATIONS 

H(I) = R (A(I)*T + B(I)*T**2 + C(I)*T**3 + D(I)*T**4) 

U(I) = R ( (A(I)-1.)*T + B(I)*T**2 + C(I)*T**3 + D(I)*T**4) 

C CP(I) = R (A(I) + 2*B(I) *T + 3*C(I) *T**2 + 4*D(I) **3) 

C 

C WHERE H(I) IS THE SENSIBLE MOLAR ENTHALPY OF SPECIES I, J/C340L 
C U(I) IS THE SENSIBLE MOLAR ENTHALPY OF SPECIES I , J/GMOL 

C CP (I) IS THE MOLAR SPECIFIC HEAT AT CONSTANT PRESSURE, J/GMOL—K 

C 

COMMON/JANAF/ A (4) ,B(4) ,C(4) ,D(4) 

DATA A/3 . 00845,2 . 94628 , 3 . 29003 , 3 . 50953/ 

DATA B/6 . 1504 IE-4 , 1 . 03 617E-3 , 2 . 694 17E-3 , 6 . 54488E-4/ 

EATA C/— i . 09800E— 7 , —3 . 39775E-7 , -8 . 66074E-7 , 9 . 41722E— 8/ 

DATA D/5. 68186E-12 , 4 . 75627E-11 , 1 . 11529E-10, -3 . 43661E-11/ 

END 
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